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1.  Introduction 


Fragment  shape  factor  7  is  defined  by  the  equation 


2/3 


=  7 


(1) 


where  Ap  is  presented  area,  m  is  mass,  and  p  is  material  density.  The  presented  (or 
projected)  area  is  the  area  that  would  be  projected  in  silhouette  on  a  screen  from  a 
distant  light  source,  and  could  very  well  change  with  orientation.  The  mass  divided 
by  the  density  is  the  fragment’s  volume.  We  raise  it  to  the  two-thirds  power  so  that 
it  has  the  same  dimensions  as  the  presented  area  (square  meters  in  SI  units).  This 
makes  the  shape  factor  7  dimensionless,  which  means  that  it  is  independent  of  any 
units.  More  importantly,  the  shape  factor  as  defined  in  Eq.  1  is  also  independent  of 
the  fragment’s  density,  so  that  the  shape  factor  of  a  steel  cube  has  the  same  value  as 
a  tungsten  cube  or  an  aluminum  cube.  It  is  also  a  specific  quantity,  independent  of 
the  volume  of  the  fragment.  This  will  prove  useful  later  when  we  compare  shape 
factors  for  different  masses. 

It  is  important  to  recognize  that  there  are  at  least  3  other  shape  factor  definitions  in 
common  use  in  the  ballistics  community: 

•  Shape  factor  s  defined  by  Ap  =  sm 2//3,  where  Ap  is  in  units  of  inch2,  m  is 
in  units  of  grains,  and  s  is  in  units  of  inch2/gr2//'3. 

•  Shape  factor  K  defined  by  Ap  =  KW 2/3,  where  Ap  is  in  units  of  ft2,  W  is 
in  units  of  pounds,  and  K  is  in  units  of  ft2  gr1//3/lb. 

•  Shape  factor  k  defined  by  Ap  =  ( m / fc)2//3,  where  Ap  is  in  units  of  inch2,  m 
is  in  units  of  grains,  and  k  is  in  units  of  gr/inch3. 


Notice  that  these  are  all  dimensional  shape  factors.  Using  the  conversion  factors  1 
inch  =  2.54  cm,  1  lb  =  7000  gr,  and  1  g  =  15.4324  gr, 


7 

•  the  conversion  between  7  and  s  is  s  =  0.025  — 


•  the  conversion  between  7  and  K  is  K  =  0.025 
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•  and  the  conversion  between  7  and  k  is  k  —  252.9—^, 

ryO/Z 

where  p  is  the  material  density  in  units  of  g/cm3.  For  example,  a  steel  cube  (p  =  7.83 
g/cm3)  with  a  random  orientation  has  a  mean  shape  factor  of 


•  7  =  3/2, 

•  s  —  0.0095  inch2/gr2//3, 

•  K  =  0.4623  ft2  gr  1//3/lb,  and 

•  k  =  1077.8  gr/inch3. 


And  a  tungsten  cube  (p  =  17.6  g/cm3)  with  a  random  orientation  has  a  mean  shape 
factor  of 


•  7  =  3/2, 

•  s  =  0.0055  inch2/gr2/3, 

•  K  =  0.2694  ft2  gr  1/3/lb,  and 

•  k  =  2422.8  gr/inch3. 


The  dimensionless  shape  factor  7  is  much  simpler  and  less  error  prone  than  the  others 
since  it  only  depends  upon  the  shape  and  orientation ,  but  is  completely  independent 
of  material  density. 

Shape  factors  for  some  common  shapes  and  orientations  can  be  worked  out  from 
the  definition  embodied  in  Eq.  1  and  simple  geometry.  Some  of  these  are  listed  in 
Table  1. 
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Table  1.  Shape  factors  of  some  common  shapes  and  orientations 


Shape 

Orientation 

Shape  Factor 

Sphere 

All 

(3/2)2/3(7t/4)1/3  «  1.209 

Cube 

Face  Forward 

1 

Edge  Forward 

1.414 

Corner  Forward 

\/3«  1.732 

Minimum 

1 

Maximum 

1.732 

Mean 

3/2 

Median 

7\/3/8  «  1.516 

3/2/1  Cuboid 

Largest  Face  Forward 

1.817 

Intermediate  Face  Forward 

0.909 

Smallest  Face  Forward 

0.606 

Minimum 

0.606 

Maximum 

2.120 

Mean 

5.5/62/3  «  1.666 

Median 

1.745 

L/D=l  Cylinder 

Face  Forward 

(tt/4)1/3  «  0.923 

Side  Forward 

(tt/4)-2/3  «  1.175 

Minimum 

(tt/4)1/3  «  0.923 

Maximum 

(7t/4)^2/3\/  (7t/4)2  +  1  ss  1.494 

Mean 

(3/2)(7r/4)1/3  «  1.384 

Median 

1.416 

Regular  Tetrahedron 

Face/Corner  Forward 

1.801 

Edge  Forward 

2.080 

Minimum 

1.471 

Maximum 

2.080 

Mean 

1.801 

Median 

1.775 

3/2/1  Ellipsoid 

Minimum 

0.732 

Mode 

1.098 

Maximum 

2.197 

Mean 

1.424 

Median 

1.368 

Shape  factor  plays  a  fundamental  role  in  ballistic  penetration.  Two  simple  examples 
will  serve  to  illustrate  this.* 

•  A  725-gr  steel  cylinder  with  a  shape  factor  7  =  0.72  ( L/D  =  1.4506),  when 
striking  a  1/4-inch  mild  steel  plate,  has  a  limit  velocity  vj,  =  1162  f/s.  If 

*A11  these  results  were  obtained  by  running  the  FATEPEN  model.12 
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we  reduce  the  shape  factor  by  50%  (7  =  0.36  and  L/ D  =  4.1029)  without 
changing  the  mass ,  then  the  limit  velocity  becomes  vl  =  943  f/s,  a  19% 
reduction  in  limit  velocity.  On  the  other  hand,  if  we  keep  the  same  shape  factor 
of  0.72,  then  the  mass  would  have  to  increase  to  1385  gr,  a  91%  mass  increase, 
in  order  to  achieve  the  same  reduction  in  limit  velocity. 

•  A  25-g  steel  cylinder  with  a  shape  factor  7  =  0.86  (L/D  =  1.11),  when 
striking  a  16-mm  face-hardened  steel  plate,  has  a  limit  velocity  vl  =  3621 
f/s.  The  same  mass  with  a  shape  factor  7  =  0.43  (L/D  =  3.14)  has  a  limit 
velocity  vl  =  3263  f/s,  a  10%  reduction  in  limit  velocity.  If  we  leave  the  shape 
factor  at  0.43,  then  the  mass  would  have  to  be  increased  to  32  g  to  get  a  limit 
velocity  of  3263  f/s,  which  represents  a  28%  increase  in  mass. 

These  examples  illustrate  that  a  decrease  in  shape  factor  is  comparable  to  an  increase 
in  striking  mass.  In  the  first  case,  a  50%  reduction  in  shape  factor  was  comparable 
to  a  91%  increase  in  mass,  and  in  the  second,  a  50%  reduction  in  shape  factor  was 
comparable  to  a  28%  increase  in  mass.  So  the  shape  factor  can  be  more  or  less 
sensitive  than  the  mass  in  influencing  the  limit  velocity.  But  the  important  point  is 
that  to  accurately  determine  penetration,  we  need  to  know  the  shape  factor  of  the 
penetrator,  much  like  the  mass.  This  is  further  illustrated  in  Fig.  1. 

1000 

800 

'c/T 

.S  600 
2 
Si 

(/> 
in 
ro 

400 

200 


1  2  3  4  5 

Shape  Factor  (-) 

Fig.  1.  Contour  plot  of  limit  velocity  as  a  function  of  shape  factor  and  mass  for  steel  fragments 
impacting  a  1/4-inch  mild  steel  plate  using  FATEPEN 
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Another  point  worth  emphasizing  is  that  penetration  depends  upon  the  instantaneous 
shape  factor  at  impact,  not  the  average  value  over  all  orientations.  Unless  the  frag¬ 
ment  is  tumbling  while  it  is  penetrating — which  is  highly  unlikely  in  metal — we 
need  to  use  the  shape  factor  at  impact.  One  may  argue  that  the  THOR  penetration 
model3  makes  use  of  the  average  presented  area,  but  a  moment’s  reflection  should 
convince  us  that  it  is  a  mistake  to  use  the  average  value.  Consider,  for  example,  a 
long  thin  cylinder.  In  face-forward  orientation,  it  is  a  very  effective  penetrator — not 
so  in  side-forward  orientation.  If  we  average  over  all  orientations,  we  may  find  that 
the  average  value  does  not  perforate.  It  would  be  a  mistake  to  conclude,  therefore, 
that  there  is  no  perforation  for  any  orientation — an  example  of  averaging  too  soon. 
Of  course,  it  is  not  necessarily  easy  to  measure  the  impact  presented  area.  For  convex 
solids,  Cauchy’s  theorem  tells  us  that  the  average  presented  area  of  a  convex  solid  is 
one-fourth  the  total  surface  area.  This  allows  us  to  compute  the  average  presented 
area  very  easily  from  the  total  surface  area  and  may  have  been  the  reason  why 
average  shape  factor  was  used  in  the  THOR  equations. 

2.  Shape  Factor  Distributions  for  5  Convex  Solids 

Now  let  us  consider  the  5  shapes  in  Table  1  that  depend  upon  orientation  (i.e., 
excluding  the  sphere),  and  let  us  consider  a  random  viewpoint  that  is  uniformly 
distributed  over  a  sphere.  We  can  imagine  the  solid  fixed  at  the  origin  while  we 
take  random  viewpoints  on  a  sphere  enclosing  the  solid,  and  from  each  viewpoint 
we  compute  the  projected  area  on  a  distant  screen  perpendicular  to  that  viewpoint. 
Analytical  formulas  for  the  projected  area  probability  density  function  (PDF)  and 
cumulative  distribution  function  (CDF)  are  known  for  each  of  these  shapes  and  are 
given  in  Appendix  A,  while  a  variety  of  methods  of  sampling  over  the  unit  sphere 
are  described  in  Appendix  B.  Recall  that  the  PDF  is  simply  the  derivative  of  the 
CDF,  and  the  CDF  always  ranges  from  0  to  1,  which  means  that  the  area  under  the 
full  range  of  the  PDF  must  equal  1. 

Some  of  the  formulas  in  Appendix  A  are  in  terms  of  the  projected  area.  To  convert 
from  the  projected  area  distribution  to  the  shape  factor  distribution,  we  make  use  of 
Eq.  1  and  the  chain  rule  to  get 

<2) 

where  /  is  the  PDF,  F  is  the  CDF,  and  V  is  the  fragment  volume.  So  we  see  that  it  is 
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easy  to  convert  a  projected  area  distribution  to  a  shape  factor  distribution  by  simple 
scaling. 

Sample  plots  of  the  PDF  and  CDF  for  each  of  the  5  shapes  are  shown  in  Subsections 
2. 1-2.5.  In  each  case  we  also  show  histograms  that  have  been  generated  with  the 
listed  simulation  code. 

2.1  Cube 

The  shape  factor  distribution  from  a  randomly  oriented  cube  can  be  simulated  with 
the  code  in  Listing  1. 


Listing  1.  cubesim.cpp 

//  cubesim.cpp:  simulate  shape  factor  distribution  of  a  cube 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

int  main(  int  argc,  char*  argv[]  )  { 

int  N  =  1000;  //  default  number  of  samples  or  specify  on  command  line 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  );  //  number  of  samples 

rng : : Random  rng; 

double  x,  y,  z,  sf; 

for  (  int  i  =  0;  i  <  N;  i++  )  { 

rng. spherical-avoidance (  x,  y,  z  ); 
sf  =  fabs(  x  )  +  fabs(  y  )  +  fabs(  z  ); 
std::cout  «  sf  «  std::endl; 

} 

return  EXIT-SUCCESS; 

} 


The  simulated  shape  factor  distribution  is  compared  to  the  analytical  formulas  in 
Fig.  2. 


Fig.  2.  Histograms  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  cube  compared  to 
analytical  formulas,  Eqs.  A-l  and  A-2  (black  curves) 
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2.2  Cuboid 


The  shape  factor  distribution  from  a  randomly  oriented  cuboid — also  known  as  a 
rectangular  parallelepiped  (RPP) — can  be  simulated  with  the  code  in  Listing  2. 


Listing  2.  rppsim.cpp 

//  rppsim.cpp 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

int  main(  int  argc,  char*  argv[]  )  { 

int  N  =  1000;  //  default  number  of  samples 

//  default  is  a  cube  with  L  =  W  =  D  =  1,  or  specify  dimensions  on  command  line 


double  L  =  1. ,  W  =  1. , 

T  =  1. ; 

if  (  argc  ==  4  )  { 

L  =  atof(  argv[l]  ) 

W  =  atof(  argv[2]  ) 

T  =  atof(  argv[3]  ) 

} 

else  if  (  argc  ==  5  )  { 

L  =  atof(  argv[l]  ) 

W  =  atof(  argv[2]  ) 

T  =  atof(  argv[3]  ) 

N  =  atoi(  argv[4]  ) 

//  number  of  samples 

} 


double  a  =  W  *  T; 
double  b  =  T  *  L; 
double  c  =  L  *  W; 
const  double  V  =  L  *  W  *  T; 

const  double  F  =  pow(  V,  -2.  /  3.  );  //  factor  to  convert  area  to  shape  factor 

rng:: Random  rng; 
double  x,  y,  z,  ap,  sf; 

for  (  int  i  =  0;  i  <  N;  i++  )  { 

rng. spherical-avoidance (  x,  y,  z  ); 

ap  =  (  a  *  fabs(  x  )  +  b  *  fabs(  y  )  +  c  *  fabs(  z  )  ); 
sf  =  ap  *  F; 

std::cout  «  sf  «  std::endl; 

} 

return  EXIT-SUCCESS; 


The  simulated  shape  factor  distribution  is  compared  to  the  analytical  formulas  in 
Fig.  3. 


Fig.  3.  Histograms  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  L  =  3,  W  —  2, 
T  —  1  cuboid  compared  to  analytical  formulas  (black  curve) 
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2.3  Cylinder 


The  shape  factor  distribution  from  a  randomly  oriented  cylinder  can  be  simulated 
with  the  code  in  Listing  3. 


Listing  3.  cylsim.cpp 


//  cylsim.cpp:  simulate  shape  factor  distribution  of  a  right  circular  cylinder 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

int  main(  int  argc,  char*  argv[]  )  { 

int  N  =  1000;  //  default  number  of  samples  or  specify  on  command  line 

double  l_d  =1.;  //  default  is  a  L/D  =  1  cylinder,  or  specify  on  command  line 

if  (  argc  ==  2  ) 

N  =  atoi(  argv[l]  );  //  number  of  samples 

else  if  (  argc  ==  3  )  { 

l_d  =  atof(  argv[l]  );  //  L/D 

N  =  atoi(  argv[2]  );  //  number  of  samples 

} 

const  double  C  =  pow(  M_PI_4  *  l_d,  -2./3.  ); 

rng:: Random  rng; 

double  th,  ph,  sf; 

for  (  int  i  =  0;  i  <  N;  i++  )  { 

rng. spherical-avoidance (  th,  ph  ); 

sf  =  C  *  (  l_d  *  sin(  th  )  +  M_PI_4  *  fabs(  cos(  th  )  )  ); 
std::cout  «  sf  «  std::endl; 

> 

return  EXIT-SUCCESS; 

> 


The  simulated  shape  factor  distribution  for  an  L/D  =  1  cylinder  is  compared  to  the 
analytical  formulas  in  Fig.  4. 


Fig.  4.  Histograms  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  L/D  —  1  cylinder 
compared  to  analytical  formulas.  Notice  the  jump  in  the  PDF  at  7  =  (7r/4)-2/3  m  1.175,  as 
predicted  (see  Appendix  A). 
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2.4  Tetrahedron 


The  shape  factor  distribution  from  a  randomly  oriented  regular  tetrahedron  can  be 
simulated  with  the  code  in  Listing  4. 


Listing  4.  tetrasim.cpp 

//  tetrasim.cpp 

#include  "Vector. h" 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

#include  <iomanip> 
using  namespace  std; 

int  main(  int  argc,  char*  argv[]  )  { 

const  double  A  =  1.  /  sqrt(  3.  ); 
const  va:: Vector  norm [4]  =  { 
va::Vector(  +A,  -A,  +A  ), 

va::Vector(  +A,  +A,  -A  ), 

va::Vector(  -A,  +A,  +A  ), 

va::Vector(  -A,  -A,  -A  ) 

}; 

rng : : Random  rng; 

double  x,  y,  z,  dotprod,  ap,  sf; 
va: : Vector  u; 

const  double  A_FACE  =  sqrt(  3.  )  /  2.; 
const  double  VOL  =  1.  /  3.; 

int  N  =  1000;  //  default  number  of  samples  or  specify  on  command  line 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  ); 

cout  «  std : : setprecision (6)  «  std:: fixed; 
for  (  int  n  =  0;  n  <  N;  n++  )  { 

ap  =  0. ; 

rng. spherical_avoidance (  x,  y,  z  ); 
u  =  va : : Vector (  x,  y,  z  ) ; 

for  (  int  i  =  0;  i  <  4;  i++  )  if  (  (  dotprod  =  u  *  norm[i]  )  >  0.  )  ap  +=  A_FACE  *  dotprod; 
sf  =  ap  *  pow(  VOL,  -2./3.  ); 
std:: cout  «  sf  «  std::endl; 

> 

return  EXIT-SUCCESS; 


The  simulated  shape  factor  distribution  for  a  regular  tetrahedron  is  compared  to  the 
analytical  formulas  in  Fig.  5. 


Fig.  5.  Histograms  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  regular  tetrahedron 
compared  to  analytical  formulas 
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2.5  Ellipsoid 


Listing  5  is  an  implementation  of  a  simulation  of  the  shape  factor  for  an  ellipsoid. 


Listing  5.  ellipsoidsim.cpp 

//  ellipsoidsim.cpp:  Simulate  the  shape  factor  a  randomly  oriented  ellipsoid 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlit» 

#include  <cmath> 


inline 

double 

min( 

double 

a,  double 

b. 

double 

c 

{ 

return 

std : 

:min( 

std : 

:min( 

a,  b  ) 

,  c  ); 

} 

inline 

double 

max( 

double 

a,  double 

b, 

double 

c 

{ 

return 

std : 

:max( 

std : 

:max( 

a,  b  ) 

,  c  ); 

> 

inline 

double 

mid( 

double 

a,  double 

b. 

double 

C 

{ 

return 

std : 

:max( 

std : 

:min( 

a,  b  ) 

,  std : 

:min(  std:: max (  a,  b  ),  c  )  ); 

} 


int  main(  int  argc,  char*  argv[]  )  { 

unsigned  int  N  =  1000;  //  default  number  of  samples  or  specify  as  4th  argument  on  command  line 

double  a  =  1. ,  b  =  1. ,  c  =  1. ;  //  default  shape  is  a  sphere 

if  (  argc  ==  4  )  {  //or  specify  the  3  dimensions  (in  any  order)  on  the  commandline 


a  =  atof(  argv[l]  ) 

b  =  atof(  argv[2]  ) 

c  =  atof(  argv[3]  ) 

} 


else  if  (  argc  ==  5  )  {  //  specify  number  of  samples  as  4th  argument 


a  =  atof(  argv[l]  ) 

b  =  atof(  argv[2]  ) 

c  =  atof(  argv[3]  ) 

N  =  atoi(  argv[4]  ) 

> 


const  double  A  =  min( 

a,  b,  c  ); 

// 

minimum  value 

const  double  B  =  mid( 

a,  b,  c  ); 

// 

intermediate  value 

const  double  C  =  max( 

a,  b,  c  ); 

// 

maximum  value 

const  double  V  =  (  4. 

/  3.  )  *  M_PI  *  A  *  B  *  C; 

// 

ellipsoid  volume 

const  double  F  =  pow( 

V,  -2.  /  3.  ); 

// 

factor  to  convert  area  to  shape  factor 

double  x,  y,  z,  X,  Y, 

Z,  ap,  sf; 

rng:: Random  rng; 

for  (  unsigned  int  n 

0;  n  <  N;  n++  )  { 

rng. spherical_avoidance (  x,  y,  z  ); 

X  =  B  *  C  *  x; 

Y  =  A  *  C  *  y; 

Z  =  A  *  B  *  z; 

ap  =  M_PI  *  sqrt(  X*X+Y  *Y+Z*  Z  ); 
sf  =  ap  *  F; 

std::cout  «  sf  «  std::endl; 

> 

return  EXIT-SUCCESS; 

> 


Running  this  code  for  a  3/2/1  ellipsoid  gives  the  results  shown  in  Fig.  6. 


Fig.  6.  Histograms  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  ellipsoid  with  a  —  1, 
b  —  2,  c  —  3.  The  solid  curve  is  a  plot  of  the  analytical  formula. 
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Running  the  code  for  a  8/2/1  ellipsoid  gives  the  results  shown  in  Fig.  7,  just  to  show 
how  the  side  length  ratios  shift  the  plots.  We  will  see  later  that  the  distribution  from 
a  range  of  ellipsoid  shapes  begins  to  resemble  the  shape  factor  distribution  from 
natural  fragments. 


Fig.  7.  Histograms  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  ellipsoid  with  a  =  1, 
b  —  2,  c  —  8.  The  solid  curve  is  the  theoretical  distribution. 


Given  the  semi-principal  axes  lengths  a,  b,  c,  of  the  ellipsoid,  where  a  <  b  <  c,  the 
projected  areas  are 


Anin  =  7T  ab,  Am  =  ttclc,  AmSLX  =  irbc,  and  V  =  -tt  abc.  (3) 

o 

Or,  if  we  know  Amin,  Am,  /lrn;ix,  where  Amin  <  Am  <  Amax,  then  the  lengths  are 


a  = 


A  ■  A 


7 TAr 


b  = 


A  A 


7 TA„ 


c  = 


A  A 

^max^m 


7 TA-r, 


and 


A  A  ■  A  A 

yr  _  ^  A  '  yrLmm^Lm^LmaX 

7T 


(4) 


Therefore,  if  we  are  given  the  minimum,  mode,  and  maximum  shape  factors,  7mim 
7m,  7max,  along  with  the  ellipsoid  volume,  then  we  can  compute 

A  ■  —  'y  ■  l/2/3  A  —'y  V21 3  A  =  'y  V2/ 3  (51 

/mm v  •>  ^*-m  /m v  •>  ^max  /max v  ?  V^/ 


and  use  Eq.  4  to  compute  the  ellipsoid  dimensions. 

Notice  carefully  the  shape  of  the  CDF  curves  for  each  of  these  shapes.  We  will  see 
that  none  of  them  resembles  the  CDF  for  a  randomly  oriented  natural  fragment, 
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which  indicates  that  there  is  no  simple  randomly  oriented  shape  that  will  suffice  as  a 
model  for  a  natural  fragment.  What  we  will  propose  instead  is  to  first  characterize 
the  probability  distribution  of  natural  fragments  and  find  a  way  to  sample  from  this 
distribution. 

3.  Shape  Factor  Distributions  for  Natural  Fragments 

Natural  fragments  resulting  from  artillery  rounds  have  been  collected  in  many 
munition  tests  over  the  years.  Unlike  the  shapes  we  have  considered  so  far,  these 
fragments  tend  to  have  highly  irregular  shapes,  an  example  of  which  is  shown  in 
Fig.  8. 


Fig.  8.  Raytraced  image  made  with  BRL-CAD  of  SF1290,  a  laser-scanned  natural  fragment 

We  make  the  assumption  that  when  a  fragment  impacts  a  target  surface,  it  has  a 
random  orientation  and  that  each  orientation  is  equally  likely.  Thus,  we  need  a  way 
to  sample  over  all  directions  in  an  unbiased  manner.  This  was  handled  in  the  past  by 
making  use  of  Platonic  solids. 

3.1  Platonic  Solids  and  Uniform  Viewing  from  All  Viewpoints 

Platonic  solids  are  3-dimensional  (3D)  regular,  convex  polyhedra,  and  there  are  only 
5,  as  shown  in  Fig.  9. 

r 

Fig.  9.  The  5  Platonic  solids,  from  left  to  right,  are  tetrahedron,  cube  or  hexahedron,  octahe¬ 
dron,  dodecahedron,  and  icosahedron 
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Regular  in  this  context  means  that  each  face  has  the  same  area,  which  is  the  key 
property  for  our  purposes.  Convex  means  that  if  we  connect  any  point  on  the  inside 
with  any  point  on  the  outside  with  a  straight  line,  then  it  will  cross  the  surface  only 
once.  The  number  of  faces  and  vertices  of  the  solids  are  listed  in  Table  2. 


Table  2.  Properties  of  the  5  Platonic  solids 


Solid 

Number  of  Faces 

Number  of  Vertices 

Tetrahedron 

4 

4 

Cube 

6 

8 

Octahedron 

8 

6 

Dodecahedron 

12 

20 

Icosahedron 

20 

12 

The  dual  of  a  Platonic  solid  is  one  in  which  the  positions  of  the  face  centers  and  the 
positions  of  vertices  are  switched — and  is  also  a  Platonic  solid.  The  tetrahedron  is 
self-dual,  and  the  hexahedron  and  octahedron  are  duals  of  one  another,  as  are  the 
dodecahedron  and  the  icosahedron.  The  vertices  of  a  Platonic  solid  are  equally  spaced 
about  a  circumscribed  sphere,  so  that  makes  them  ideal  candidates  for  unbiased 
projected  area  viewpoints. 

The  Icosahedron  Gage 4-6  is  a  measuring  instrument  that  uses  as  viewpoints  the 
vertices  of  both  the  icosahedron  and  the  dodecahedron.  This  gives  it  12  +  20  =  32 
viewpoints,  but  for  projected  areas,  half  of  these  viewpoints  are  redundant,  so  that 
leaves  16  viewpoints.*  Table  3  lists  the  angles  of  the  Icosahedron  Gage. 


*Notice,  however,  that  when  we  do  this,  we  no  longer  have  equal  spacing  between  viewpoints. 
When  a  Platonic  solid  and  its  dual  are  combined,  some  points  have  3  closest  neighbor  vertices  while 
other  points  have  5.  So  the  rotational  symmetry  is  spoiled.  There  is  simply  no  way  to  distribute  more 
than  20  points  over  the  unit  sphere  and  maintain  equal  spacing  between  nearest  neighbors.  If  it  were 
possible,  then  we  would  have  a  new  Platonic  solid,  but  it  has  been  proven  that  there  are  only  5. 
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Table  3.  Sequence  of  viewing  angles  in  Icosahedron  Gage 


Position 

Platonic  Solid 

Azimuthal  Angle 

(°) 

Elevation  Angle 

(°) 

1 

Icosahedron 

0 

90 

2 

Dodecahedron 

0 

52.6226 

3 

72 

52.6226 

4 

144 

52.6226 

5 

216 

52.6226 

6 

288 

52.6226 

7 

Icosahedron 

324 

26.5651 

8 

36 

26.5651 

9 

108 

26.5651 

10 

180 

26.5651 

11 

252 

26.5651 

12 

Dodecahedron 

288 

10.8123 

13 

0 

10.8123 

14 

72 

10.8123 

15 

144 

10.8123 

16 

216 

10.8123 

Projected  area  measurements  have  been  performed  with  early  versions  of  the  Icosa¬ 
hedron  Gage  since  the  1940s.  The  instrument  that  is  used  today  is  coupled  to  a 
personal  computer,  which  greatly  automates  the  process.6,7 

3.2  Natural  Fragments  from  Artillery  Rounds 

Close  to  900  steel  fragments  from  artillery  rounds  have  been  collected.8-11  Each 
fragment  mass  was  measured  along  with  16  projected  areas  with  an  Icosahedron 
Gage,  which  allows  us  to  compute  16  shape  factor  values  for  each  fragment. 

If  we  plot  the  mean  shape  factor  (averaged  over  the  16  individual  measurements)  as 
a  function  of  fragment  mass,  we  find  that  there  is  essentially  no  correlation  between 
fragment  mass  and  average  shape  factor,  as  shown  in  Fig.  10. 
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Fig.  10.  Fragments  from  122-mm,  152-mm,  and  155-mm  artillery  rounds.  Mean  shape  factors 
are  computed  from  the  16  viewpoints  of  the  Icosahedron  Gage.  The  dashed  line  is  a  least- 
squares  fit  to  the  data  with  a  slope  of  very  nearly  zero  (0.0018)  indicating  essentially  no  corre¬ 
lation  between  average  shape  factor  and  mass. 

Since  there  is  essentially  no  correlation  between  shape  factor  and  mass,*  we  are 
justified  in  pooling  all  of  the  shape  factors,  which  then  gives  us  a  sample  size  of 
898  x  16  =  14,  368  shape  factors. 

So,  although  we  started  out  with  only  16  shape  factors  for  each  irregular  fragment, 
we  are  able  to  exploit  the  fact  that  shape  factor  is  independent  of  mass  to  effectively 
come  up  with  over  14,000  shape  factor  measurements.  We  find  that  the  resulting 
shape  factor  distribution  closely  approximates  a  lognormal  distribution,  as  shown  in 
Fig.  11. 


Fig.  11.  Histogram  of  shape  factor  PDF  and  CDF  for  898  artillery  fragments.  The  black  curve 
is  the  maximum  likelihood  fit  to  the  lognormal  distribution  with  //  =  0.597  and  a  —  0.341. 

*For  example,  we  might  have  expected  that  smaller  fragments  would  be  more  compact  than 
larger  fragments,  but  that  is  not  what  we  see  in  the  data. 
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(6) 


(7) 


The  maximum  likelihood  estimation*  of  parameters  gives  /i  =  0.596514  and  a  = 
0.340874.  The  geometric  mean  is  =  eM  =  1.81578  and  the  geometric  standard 
deviation  is  ag  =  ea  =  1.40618.  '  The  comparison  between  the  lognormal  fit  and  the 
data  is  summarized  in  Table  4. 


Table  4.  Comparison  of  lognormal  fit  to  artillery  data 


Whereas  a  normal  distribution  has  the  property  that  p  ±  x  are  equally  likely,  a 

2  1  2 

lognormal  distribution  has  the  property  that  xe^~a  and  are  equally  likely, 

for  any  value  x  ^  0.  See  Appendix  C  for  some  more  properties  of  the  lognormal 
distribution. 

3.3  Natural  Fragments  from  Spall 

Spall  fragments*  were  also  collected  in  Celotex  and  subsequently  measured  for 
mass  and  projected  area  with  an  Icosahedron  Gage.  In  this  case  there  were  250  spall 
fragments,  giving  250  x  16  =  4000  shape  factors.  These  data  are  displayed  in  Fig.  12 
and  again  compared  to  a  lognormal  fit. 


*  Maximum  likelihood  estimation  is  easy  to  perform  for  a  lognormal  distribution  since  //  and  a 
are  respectively  the  mean  and  the  standard  deviation  of  the  logs  of  the  data. 

^The  geometric  mean  and  geometric  standard  deviation  make  it  easy  to  summarize  the  distribution. 
Thus,  68%  is  contained  in  \'yg(Jg1, 7ga9]  and  95%  in  [7s(j“  2,  7scr2]. 

Spall  fragments  are  pieces  of  armor  that  are  broken  off  during  penetrator  impact. 
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Fig.  12.  Histogram  of  shape  factor  PDF  and  CDF  for  spall  fragments.  The  black  curve  is  a  fit 
to  the  lognormal  distribution. 


Maximum  likelihood  estimation  gives  /i  =  0.456095  and  a  =  0.479386.  The 
geometric  mean  is  7 9  =  1.5779  and  the  geometric  standard  deviation  is 

ag  =  ea  =  1.61508.  The  comparison  between  the  lognormal  fit  and  the  data  is 
summarized  in  Table  5. 


Table  5.  Comparison  of  lognormal  fit  to  spall  data 


Statistic 

Lognormal  Fit 

Spall  Data 

Median 

=  1.58 

1.50 

Mean 

efi+a2/2  =  177 

1.79 

Mode 

e v-*2  =  1.25 

1.25 

If  we  look  back  and  compare  the  PDF  for  the  simple  shapes  of  cube,  cuboid,  cylinder, 
tetrahedron,  and  ellipsoid  (Figs.  1-6)  to  the  lognormal,  it  is  clear  that  none  of  these 
shapes  come  close.  Nor  are  we  likely  to  find  any  simple  shape  that  will  reproduce 
the  distribution  of  shape  factors  from  natural  fragments  by  merely  randomizing  the 
orientation.  However,  we  have  found  that  the  lognormal  distribution  offers  a  lot  of 
promise  of  simulating  the  shape  factor.  It  is  not  necessary  that  we  have  a  specific 
shape  in  order  to  simulate  the  shape  factor;  we  only  need  to  sample  a  lognormal 
distribution  each  time  we  need  a  sample  of  the  shape  factor  (for  example,  at  impact 
with  a  target). 

Listing  6  is  a  simple  program  that  will  generate  the  lognormally  distributed  shape 
factors. 
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Listing  6.  lognormal.cpp 


//  lognormal.cpp 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

#include  <chrono> 

#include  <random> 

int  main(  void  )  { 

//  default  values  for  the  shape  factor  lognormal  distribution  from  122mm,  152mm  and  155mm  artillery 
double  mu  =  0.590494;  //  these  two  parameters  characterize  the  lognormal  shape  factor  distribution 

double  sigma  =  0.323433;  //  with  mode  =  1.63,  median  =  1.80  and  mean  =  1.90 

unsigned  int  seed  =  std :: chrono :: high_resolution_clock: : now( ). time_since_epoch( ). count () ; 
std::mtl9937  rng(  seed  );  //  Mersenne  Twister  engine 

std : : lognormal_distribution<double>  lognormal!  mu,  sigma  );  //  lognormal  shape  factor  distribution 

const  int  N  =  10000; 

for  (  int  i  =  0;  i  <  N;  i++  )  std::cout  «  lognormal!  rng  )  «  std::endl; 
return  EXIT-SUCCESS; 

} 


3.4  Shape  Factor  Computation  from  Laser-Scanned  Fragments 

Another  technique  that  has  been  used  more  recently  for  measuring  fragment  shape 
is  laser  scanning.1  This  will  generate  a  facetized  solid  in  stereolithography  (STL) 
format.12  An  example  is  shown  in  Fig.  13. 


Fig.  13.  Laser-scanned  natural  fragment  showing  mesh  of  470,988  triangles  covering  the  sur¬ 
face.  It  is  obvious  that  the  fragment  is  not  convex,  which  means  that  the  many  hidden  surfaces 
need  to  be  accounted  for  when  computing  the  projected  area. 


Listing  7  shows  the  format  for  an  STL  file. 

Listing  7.  stlformat.cpp 

//  stlformat.cpp 
solid  name 

facet  normal  nl  n2  n3 
outer  loop 

vertex  vlx  vly  viz 
vertex  v2x  v2y  v2z 
vertex  v3x  v3y  v3z 
endloop 
endfacet 
endsolid 
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Thus,  each  facet  is  a  triangle,  specified  by  4  vectors: 


•  an  outward  normal  vector,  which  follows  the  right-hand  rule*  and 

•  one  vector  for  each  of  its  3  vertices 


For  example,  Listing  8  is  the  STL  file  (in  ASCII  format^)  for  a  regular  tetrahedron. 


Listing  8.  tetrahedron.stla 


solid  TETRAHEDRON 

facet  normal  0.57735  -0.57735  0.57735 
outer  loop 

vertex  001 
vertex  100 
vertex  111 
endloop 
endfacet 

facet  normal  0.57735  0.57735  -0.57735 
outer  loop 

vertex  111 
vertex  100 
vertex  010 
endloop 
endfacet 

facet  normal  -0.57735  0.57735  0.57735 
outer  loop 

vertex  111 
vertex  010 
vertex  001 
endloop 
endfacet 

facet  normal  -0.57735  -0.57735  -0.57735 
outer  loop 

vertex  001 
vertex  010 
vertex  100 
endloop 
endfacet 

endsolid  TETRAHEDRON 


The  format  is  somewhat  redundant  in  that  the  outward  normal  can  be  computed 
from  the  3  vertices  and  is  therefore  not  strictly  required.  Indeed,  some  software 
applications  do  not  specify  the  outward  normal  in  the  STL  file,  so  one  must  be 
careful  to  check  the  normals.  Listing  9  is  a  program  to  read  in  an  STL  file  in  ASCII 
format  and  print  summary  information. 

Listing  9.  stl.a.cpp 

//  stl.a.cpp:  reads  in  an  ASCII  STL  decription  of  a  solid  and  computes  volume  and  surface  area 

#include  <iostream> 

#include  <cstdlib> 

#include  <string> 

#include  <cmath> 

#include  <iomanip> 


*The  right-hand  rule  specifies  that  if  the  fingers  of  the  right  hand  curl  in  the  direction  of  the 
corners  of  the  triangle  from  its  first  to  second  to  third  point,  then  the  thumb  will  be  pointing  in  the 
direction  of  the  outward  normal. 

There  is  also  a  binary  STL  format,  which  we  shall  use  in  Listing  10. 
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inline  double  square(  double  x  )  {  return  x  *  x;  } 


int  main(  void  )  { 

std::string  solid,  name,  facet,  normal,  outer,  loop,  vertex,  endloop,  endfacet,  endsolid; 
double  nx,  ny,  nz,  vlx,  vly,  viz,  v2x,  v2y,  v2z,  v3x,  v3y,  v3z; 
double  ax,  ay,  az,  bx,  by,  bz; 

std::cin  »  solid  »  name; 

int  n_facets  =  0; 

double  total_area  =  0.,  volume  =  0.,  mean_area; 

while  (  std::cin  »  facet  »  normal  »  nx  »  ny  »  nz  )  {  //  for  each  facet 

std::cin  »  outer  »  loop; 

std::cin  »  vertex  »  vlx  »  vly  »  viz; 

std::cin  »  vertex  »  v2x  »  v2y  »  v2z; 

std::cin  »  vertex  »  v3x  »  v3y  »  v3z; 

std::cin  »  endloop; 
std::cin  »  endfacet; 

n_facets++; 

//  compute  twice  the  area  of  the  triangle  specified  by  its  three  vertices 

ax  =  v2x  -  vlx; 

ay  =  v2y  -  vly; 

az  =  v2z  -  viz; 

bx  =  v3x  -  vlx; 

by  =  v3y  -  vly; 

bz  =  v3z  -  viz; 

total_area  +=  sqrt(  square(  ay  *  bz  -  az  *  by  )  +  square(  az  *  bx  -  ax  *  bz  )  +  square(  ax  *  by  -  ay  *  bx  )  ); 

//  compute  six  times  the  volume  of  the  tetrahedron  formed  by  the  given  triangle  and  the  origin 
//  note  that  this  is  an  oriented  volume,  which  could  be  positive  or  negative, 

//  and  the  origin  is  completely  arbitrary  so  might  as  well  use  Vector(0,0 ,0) 

volume  +=  vlx  *  (  v2y  *  v3z  -  v2z  *  v3y  )  +  vly  *  (  v2z  *  v3x  -  v2x  *  v3z  )  +  viz  *  (  v2x  *  v3y  -  v2y  *  v3x  ) ; 

} 

total-area  /=  2.; 
volume  /=  6 . ; 

mean_area  =  total-area  /  4.;  //  from  Cauchy's  theorem 


std : 

:  cout 

« 

std : : setprecision (6)  « 

stc 

: : fixed ; 

std : 

:  cout 

« 

"Solid  name  =  " 

« 

name  «  std : : endl; 

std : 

:  cout 

« 

"Number  of  facets  =  " 

« 

n_facets  «  std:: endl; 

std : 

:  cout 

« 

"Total  surface  area  =  " 

« 

total_area  «  std:: endl; 

std : 

:  cout 

« 

"Volume  =  " 

« 

volume  «  std:: endl; 

std : 

:  cout 

« 

"Mean  surface  area  =  " 

« 

total_area  /  4.  «  "  (from  Cauchy's 

theorem)"  «  std:: endl; 

std : 

:  cout 

« 

"Mean  shape  factor  =  " 

« 

mean_area  *  pow(  volume,  -2./3.  )  « 

"  (from  Cauchy's  theorem)"  «  std:: endl; 

return  EXIT-SUCCESS; 


3.4.1  Area  Contribution  from  Each  Facet 

Implicit  in  the  specification  of  the  triangle  is  an  origin  for  its  vector  vertices,  but  the 
origin  itself  is  not  explicitly  specified  in  the  STL  file.  It  is  here  that  the  orientation 
of  the  triangle  comes  to  our  aid,  since  it  turns  out  that  the  origin  is  not  necessary  to 
compute  the  triangle’s  area. 

Using  the  fact  that  the  area  of  a  triangle  is  one-half  the  area  of  the  corresponding 
parallelogram,  the  area  of  the  triangle  can  be  written  as 


where  n  is  the  outward  normal,  and  v1;  v2,  and  v3  are  the  3  vector  vertices.  Expand- 
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ing  this  out  provides  another  formula: 

A  =  ^  n  •  [vi  x  v2  +  v2  x  v3  +  v3  x  vi],  (9) 

Interpreting  this  geometrically,  it  is  easy  to  see  that  each  of  the  3  terms  is  positive  if 
the  projection  of  the  origin  lies  inside  the  triangle,  but  that  one  of  these  terms  will  be 
negative  if  the  projected  origin  lies  outside  the  triangle.  The  net  result  is  always  the 
triangle  area.  Equation  8  is  slightly  more  efficient  in  computer  code  than  Eq.  9  since 
it  requires  2  vector  subtractions  and  only  1  vector  cross  product,  as  opposed  to  2 
additions  and  3  cross  products.  Since  the  vector  cross  product  in  Eq.  8  points  in  the 
same  direction  as  the  outward  normal,  the  area  of  the  triangle  can  also  be  computed 
from  the  magnitude  of  the  cross  product: 

A  =  ^  ||(v2  -  vi)  x  (v3  -  vi) 1 1 .  (10) 

We  have  implemented  3D  vectors  and  their  associated  algebra  directly  into  software 
as  a  C++  class,  so  that  Eqs.  8  and  10  can  be  coded  directly.  Alternatively,  these 
formulas  can  be  reduced  to  scalar  formulas  by  using 

'i  j  k 

a  x  b  =  det  ax  ay  az 
bx  by  bz 

=  \{aybz  -  azby )  +  j (azbx  -  axbz )  +  k (axby  -  aybx ),  (11) 

where  i,  j,  and  k  are  3  unit  vectors  along  the  x,  y,  and  z  axes,  respectively.  Therefore, 
the  area  can  be  written  as 


or,  since  n  is  a  unit  vector  normal  to  the  facet, 

1 

A  -^\fix[®jybz  azby)  T  tty(cizbx  ttxbz^j  ~t~  tiz(^axby  ciybx^.  (13) 

3.4.2  Volume  Contribution  from  Each  Facet 

Each  facet,  combined  with  another  fixed  point  d,  makes  up  a  tetrahedron.  Let 
a  =  vi  —  d,  b  =  v2  —  d,  and  c  =  v3  —  d,  then  the  volume  of  this  tetrahedron  is  the 
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scalar  or  triple  product 


V=-a-(bxc).  (14) 

o 

The  fixed  point  d  is  arbitrary,  so  we  might  as  well  choose  the  origin,  in  which  case 
the  formula  for  the  volume  is  simply 


V  =  -  vi  •  (v2  x  v3). 
o 


(15) 


This  triple  product  is  symmetric  in  all  3  vertices,  but  the  order  is  important  since 
this  is  an  oriented  volume.  Once  again,  we  can  code  this  directly  in  the  C++  Vector 
Class,  or  we  can  make  use  of  the  following  formula: 


a  x  b 


det 


Qjx  Ciy  QjZ 

bx  by  bz 


(16) 


so  that 

V  =  ^\vix(v2yv3z  -  v2zv3y)  +  vly(v2zv3x  -  v2xv3z)  +  vlz{y2xv3y  -  v2yv3x)\.  (17) 


We  purposely  did  not  take  the  absolute  value  of  this  expression  because  these  are 
oriented  volumes.  We  do  not  really  know  the  origin  for  the  laser  scanner.  It  could  be 
inside  or  outside  the  fragment.  If  it  is  outside  the  fragment,  then  some  of  the  facets 
will  contribute  a  negative  volume  to  the  total  volume.  So,  although  the  individual 
volume  contributions  depend  upon  the  choice  of  the  vector  d,  the  total  volume  does 
not. 


3.4.3  Total  Surface  Area  and  Total  Volume 

By  adding  up  individual  contributions  to  the  surface  area  and  volume  from  Eqs.  8 
and  15,  we  can  compute  the  total  surface  area  and  total  volume  of  the  fragment. 
These  formulas  hold  whether  or  not  the  solid  is  convex. 


3.4.4  Projected  Area 

Let  the  viewing  direction  be  specified  by  a  unit  vector  u.  Then  the  area  of  the  triangle 
projected  on  a  plane  that  is  perpendicular  to  u  is  given  by 

Ap  =  ^  u  •  [(v2  —  v3)  x  (v3  —  v3)]  if  and  only  if  u  •  n  >  0.  (18) 
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The  requirement  that  u  •  n  >  0  is  necessary  so  that  we  only  compute  the  projected 
area  of  one  side  of  the  triangle,  not  both  sides.  If  the  solid  is  convex,  then  the  total 
projected  area  from  the  given  viewpoint  u  is  the  simple  summation  from  all  the 
triangles.  These  fragments  are  not  convex,  however,  which  means  that  they  have 
hidden  surfaces.  Not  taking  this  into  account  will  overestimate  the  presented  area  of 
the  fragment.  To  properly  account  for  hidden  surfaces,  we  can  use  raytracing  from 
BRL-CAD  to  compute  the  presented  area  from  a  given  direction.  It  is  useful  to  know 
how  much  the  fragments  deviate  from  convexity,  so  we  also  compute  the  presented 
area  by  simply  adding  the  contributions  according  to  Eq.  18. 


The  computation  of  the  total  projected  area  obtained  by  summing  the  contributions 
from  each  triangle,  using  Eq.  18,  is  fast  enough  that  we  have  used  10,000  points  over 
the  unit  sphere.*  But  we  also  know  that  these  fragments  are  not  convex,  so  we  have 
also  used  1000  points  distributed  over  the  unit  sphere  to  raytrace  a  projected  area 
that  fully  accounts  for  hidden  surfaces.  Once  the  projected  area  has  been  computed, 
the  (dimensionless)  shape  factor,  7,  is  then  computed  from  7  =  APV~2^. 


As  an  approximation  to  the  true  projected  area,  we  can  compute  the  projected  area 
of  the  STL  solid  without  accounting  for  hidden  surfaces  and  the  fact  that  the  solid  is 
not  convex,  as  in  Listing  10. 


Listing  10.  stl.b.cpp 

//  stl.b.cpp:  reads  in  a  binary  STL  decription  of  a  solid  and  computes  volume  and  surface  area 

#include  "Vector. h" 

#include  "Random. h" 

#include  <fstream> 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

#include  <iomanip> 

#include  <vector> 

inline  double  square(  double  x  )  {  return  x  *  x;  } 

struct  facet  { 

float  nx,  ny,  nz; 
float  vlx,  vly,  viz; 
float  v2x,  v2y,  v2z; 
float  v3x,  v3y,  v3z; 
unsigned  int  byte_count; 

}; 

int  main(  int  argc,  char*  argv[]  )  { 

char  *  file; 
if  (  argc  ==  2  )  { 
file  =  argvfl] ; 

} 

else  { 

std::cerr  «  "Usage:  "  «  argv[0]  «  "  stl_binary_inputfile  >  outputfile"  «  std::endl; 
exit (  1  ) ; 

} 

char  header [100]  = 


*See  Appendix  B  for  sampling  strategies. 
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34  unsigned  long  int  n_facets  =  0; 

35  facet  tri; 

36 

37  std : : if stream  fin; 

38  fin. open (  file,  std::ios::in  |  std :: ios :: binary  ); 

39  if  (  fin . bad ( )  )  { 

40  std::cerr  «  "Error  in  opening  file  "  «  file  «  std::endl; 

41  exit (  1  ) ; 

42  } 

43  std::cerr  «  "Contents  of  "  «  file  «  "  «  std::endl; 

44  fin. read (  (  char  *  )  header,  80  ); 

45  fin. read (  (  char  *  )  &n_facets ,  4  ); 

46  std::cerr  «  "Header  =  "  «  header  «  std::endl; 

47  std::cerr  «  "There  are  "  «  n_facets  «  "  triangles  in  "  «  file  «  "  «  std::endl; 

48 

49  double  nx,  ny,  nz,  vlx,  vly,  viz,  v2x,  v2y,  v2z,  v3x,  v3y,  v3z; 

50  double  ax,  ay,  az,  bx,  by,  bz; 

51 

52  double  total_area  =  0.,  volume  =  0.,  mean_area,  a_facet; 

53  std::vector<  va:: Vector  >  norm; 

54  std::vector<  double  >  area; 

55  for  (  unsigned  int  i  =  0;  i  <  n_facets;  i++  )  {  //  for  each  facet 

56 

57  fin. read (  (  char  *  )  &tri,  50  ); 

58 

59  //  note  that  we  don't  need  the  normal  vector  to  compute  the  area  or  volume 


60 

nx  = 

double( 

tri. 

nx  ); 

61 

ny  = 

double( 

tri. 

ny  ); 

62 

nz  = 

double( 

tri. 

nz  ); 

63 

64 

norm 

.push_back(  va::Vector(  nx,  ny,  nz  )  );; 

65 

66 

vlx 

=  double ( 

tri 

.vlx  ) 

67 

vly 

=  doublet 

tri 

•  vly  ) 

68 

viz 

=  doublet 

tri 

.viz  ) 

69 

v2x 

=  double( 

tri 

v2x  ) 

70 

v2y 

=  doublet 

tri 

v2y  ) 

71 

v2z 

=  doublet 

tri 

v2z  ) 

72 

v3x 

=  doublet 

tri 

v3x  ) 

73 

v3y 

=  doublet 

tri 

v3y  ) 

74 

v3z 

=  doublet 

tri 

v3z  ) 

75 

76  //  compute  the  area  of  the  triangle  specified  by  its  three  vertices 

77  ax  =  v2x  -  vlx; 

78  ay  =  v2y  -  vly; 

79  az  =  v2z  -  viz; 

80  bx  =  v3x  -  vlx; 

81  by  =  v3y  -  vly; 

82  bz  =  v3z  -  viz; 

83  a_facet  =  0.5  *  sqrt(  square(  ay  *  bz  -  az  *  by  )  +  square(  az  *  bx  -  ax  *  bz  )  +  square(  ax  *  by  -  ay  *  bx  )  ); 

84  total-area  +=  a_facet; 

85 

86  area . push_back(  a_facet  ); 

87 

88  //  compute  six  times  the  volume  of  the  tetrahedron  formed  by  the  given  triangle  and  the  origin 

89  //  note  that  this  is  an  oriented  volume,  which  could  be  positive  or  negative, 

90  //  and  the  origin  is  completely  arbitrary  so  might  as  well  use  Vector(O,0,0) 

91  volume  +=  vlx  *  (  v2y  *  v3z  -  v2z  *  v3y  )  +  vly  *  (  v2z  *  v3x  -  v2x  *  v3z  )  +  viz  *  (  v2x  *  v3y  -  v2y  *  v3x  ) ; 

92  > 

93  volume  /=  6.; 

94  mean_area  =  total_area  /  4.; 

95  const  double  F  =  pow(  volume,  -2./3.  ); 

96 

97  std:: clog  «  std : : setprecision (6)  «  std:: fixed; 

98  std:: clog  «  "Number  of  facets  =  "  «  n_facets  «  std::endl; 

99  std:: clog  «  "Total  surface  area  =  "  «  total_area  «  std::endl; 

100  std:: clog  «  "Volume  =  "  «  volume  «  std::endl; 

101  std::clog  «  "Mean  surface  area  =  "  «  total_area  /  4.  «  "  (from  Cauchy's  theorem)"  «  std::endl; 

102  std::clog  «  "Mean  shape  factor  =  "  «  mean_area  *  F  «  "  (from  Cauchy's  theorem)"  «  std::endl; 

103 

104  const  int  N  =  100000; 

105  rng : : Random  rng ; 

106  double  dotprod,  sf,  ap,  x,  y,  z; 

107  va: : Vector  u; 

108  std::cout  «  std :: setprecision (6)  «  std:: fixed; 

109  for  (  int  n  =  0;  n  <  N;  n++  )  { 

110 

111  ap  =  0.; 

112  rng. spherical_avoidance(  x,  y,  z  ); 

113  u  =  va::Vector(  x,  y,  z  ); 

114 

115  for  (  unsigned  int  i  =  0;  i  <  n_facets;  i++  )  if  (  (  dotprod  =  u  *  norm[i]  )  >  0  )  ap  +=  area[i]  *  dotprod; 

116  sf  =  ap  *  F; 

117  std::cout  «  sf  «  std::endl; 

118  > 
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119 

120  return  EXIT_ SUCCESS; 

121  } 


When  this  is  applied  to  the  fragment  (SF1290)  shown  in  Fig.  8,  we  get  the  histogram 
of  shape  factors  displayed  in  Fig.  14. 


Fig.  14.  Shape  factor  histograms  generated  by  STL  description  of  SF1290.  No  account  has 
been  taken  for  hidden  surfaces  so  this  is  an  overestimate  of  the  shape  factor. 


Notice  the  similarity  to  the  shape  factor  distribution  from  an  ellipsoid  (cf.  Figs.  6 
and  7).  Now  notice  what  happens  when  we  sample  from  all  the  15  fragments  that 
were  scanned7  (Fig.  15). 


Fig.  15.  Shape  factor  histograms  from  15  fragments  in  STL  form.  The  routine  used  to  gen¬ 
erate  this,  stl.b.cpp,  does  not  account  for  hidden  surfaces,  and  since  the  fragments  are  not 
convex,  this  is  an  overestimate.  The  black  curve  is  the  lognormal  fit.  The  maximum  likelihood 
estimate  is  fi  —  0.607588  and  cr  —  0.275688,  which  gives  a  median  eJ1  —  1.71746,  mean 
g/i+o-2/2  _  1.7949,  and  mode  e^-0-2  =  1.57247. 


Once  again  we  see  that  the  resulting  distribution  begins  to  approximate  a  lognormal 
distribution. 
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The  next  step  is  to  take  account  of  the  hidden  surfaces  to  come  up  with  a  more 
accurate  computation  of  the  true  projected  area.  The  computer  code  in  Listing  11 
was  used  for  this  purpose,  and  we  get  the  results  shown  in  Fig.  16. 


Fig.  16.  Shape  factor  histograms  for  15  fragments  using  stlarea.  This  computer  code  com¬ 
putes  the  projected  area  and  does  account  for  hidden  surfaces,  much  like  the  Icosahe¬ 
dron  Gage.  The  black  curve  is  the  lognormal  fit.  The  maximum  likelihood  estimate  is 
/x  =  0.537  and  <x  =  0.297,  which  gives  a  median  =  1.71,  mean  e^+CT  /2  =  1.79,  and 
mode  =  1.57. 


Listing  11.  stlarea.cpp 

//  stlarea.cpp:  Read  an  STL  binary  file  and  compute  presented  areas  and  shape  factors 
//  using  viewpoints  from  the  spiral  distribution  over  the  unit  sphere. 

//  This  implementation  computes  a  true  presented  area,  similar  to  raytracing, 

//  by  imposing  a  grid  over  the  bounding  box  from  each  viewpoint  and  then 

//  checking  if  the  grid  point  lies  in  at  least  one  of  the  projected  facets. 

//  Note:  Little  endian  is  assumed  in  STL  so  may  need  to  convert  when  porting  to  another  computer. 
//  R.  Saucier,  September  2011 

#include  "Rotation. h" 

#include  <fstream> 

#include  <iostream> 

#include  <iomanip> 

#include  <cstdlib> 

#include  <string> 

#include  <vector> 

//using  namespace  std; 

struct  facet  { 

float  nx,  ny,  nz; 
float  vlx,  vly,  viz; 
float  v2x,  v2y,  v2z; 
float  v3x,  v3y,  v3z; 
unsigned  int  byte_count; 

}; 


va::Vector  normal!  const  va::Vector&  vl,  const  va::Vector&  v2,  const  va::Vector&  v3  )  {  //  returns  outward  normal  of 

oriented  triangle 

return  unit(  (  v2  -  vl  )  A  (  v3  -  vl  )  ); 

} 

//  returns  twice  the  area  of  the  triangle  specified  by  its  three  vertices 

double  area2(  const  va:: Vectors  vl,  const  va:: Vectors  v2,  const  va:: Vectors  v3  )  {  //  one  way  to  compute  the  area 

return  mag(  (  v2  -  vl  )  *  (  v3  -  vl  )  ); 

} 

//  returns  twice  the  area  of  the  triangle  specified  by  its  normal  vector  and  its  three  vertices 

double  area2(  const  va::VectorS  n,  const  va::VectorS  vl,  const  va::VectorS  v2,  const  va::VectorS  v3  )  {  //  alternative 

way  to  compute  the  area 

return  n  *  (  (  v2  -  vl  )  A  (  v3  -  vl  )  ); 

} 
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43 

44 

//  returns  six  times  the  volume  of  the  tetrahedron  formed  by  the  given  triangle  and  the  origin 

45 

//  note  that  this  is  an  oriented  volume,  which  could  be  positive  or  negative, 

46 

//  and  the  origin  is  completely  arbitrary  so  might  as  well  use  Vector(0,0,0) 

47 

double  volume6(  const  va::Vector&  vl,  const  va::Vector&  v2,  const  va::Vector&  v3  )  { 

//  returns  the  oriented  volume 

48 

49 

return  vl  *  (  v2  /v  v3  ) ; 

50 

> 

51 

52 

inline  double  min(  double  a,  double  b,  double  c  )  {  return  std::min(  std::min(  a,  b  ) 

c  );  } 

53 

54 

inline  double  max(  double  a,  double  b,  double  c  )  {  return  std::max(  std::max(  a,  b  ) 

c  );  } 

55 

56 

//  fast  point-in-triangle  test  (Ref:  www.blackpawn.com/texts/pointinpoly/default.html 

57 

bool  inside (  const  va::Vector&  p,  const  va:: Vectors  a,  const  va:: Vectors  b,  const  va: 

Vector&  c  )  { 

58 

59 

//  compute  vectors 

60 

va: : Vector  v0  =  c  -  a; 

61 

va: : Vector  vl  =  b  -  a; 

62 

va: : Vector  v2  =  p  -  a; 

63 

64 

//  compute  dot  products 

65 

double  dot00  =  v0  *  v0; 

66 

double  dot01  =  v0  *  vl; 

67 

double  dot02  =  v0  *  v2; 

68 

double  dotll  =  vl  *  vl; 

69 

double  dotl2  =  vl  *  v2; 

70 

71 

//  Compute  barycentric  coordinates 

72 

double  w  =  1.  /  (  dot00  *  dotll  -  dot01  *  dot01  ); 

73 

double  u  =  (  dotll  *  dot02  -  dot01  *  dotl2  )  *  w; 

74 

double  v  =  (  dot00  *  dotl2  -  dot01  *  dot02  )  *  w; 

75 

76 

//  Check  if  point  is  in  triangle 

77 

return  (u>0)&&(v>0)&&(u+v<l); 

78 

> 

79 

80 

double  area_projected (  const  std : : vector<va : : Vector>S  facet_vl,  const  std : : vectorcva : 

Vector>S  facet_v2,  const  std::vecto 

<va : : Vector>&  facet_v3,  const  va:: Vectors  u,  int  n_dim  )  { 

81 

82 

const  int  N-FACETS  =  facet_vl . size( ) ; 

83 

static  const  va::Vector  I(  1.,  0.,  0.  ),  J(  0.,  1.,  0.  ),  K(  0.,  0.,  1.  ); 

84 

va:: Vector  vl,  v2,  v3; 

85 

double  x_min  =  l.e36,  x_max  =  -l.e36,  y_min  =  l.e36,  y_max  =  -l.e36; 

86 

double  vlx,  v2x,  v3x,  vly,  v2y,  v3y; 

87 

double  xmin,  xmax,  ymin,  ymax; 

88 

89 

va: : Vector  i  =  I; 

90 

va: : Vector  j  =  J; 

91 

va:: Rotation  R; 

92 

va::Vector  minus_u  =  -1.  *  u; 

93 

if  (  u  !=  K  &&  minus_u  !=  K  )  { 

94 

R  =  va :: Rotation (  K,  u  ); 

95 

i  =  R  *  I; 

96 

j  =  R  *  J ; 

97 

> 

98 

99 

//  compute  the  bounding  box  in  the  plane  perpendicular  to  u 

100 

for  (  int  n  =  0;  n  <  N_FACETS;  n++  )  { 

101 

102 

vl  =  facet_vl[n] ; 

103 

v2  =  facet_v2[n]; 

104 

v3  =  facet_v3[n]; 

105 

vlx  =  vl  *  i; 

106 

v2x  =  v2  *  i; 

107 

v3x  =  v3  *  i; 

108 

vly  =  vl  *  j ; 

109 

v2y  =  v2  *  j ; 

110 

v3y  =  v3  *  j ; 

111 

xmin  =  min(  vlx,  v2x,  v3x  ); 

112 

xmax  =  max(  vlx,  v2x,  v3x  ); 

113 

ymin  =  min(  vly,  v2y,  v3y  ); 

114 

ymax  =  max(  vly,  v2y,  v3y  ); 

115 

x_min  =  xmin  <  x_min  ?  xmin  :  x_min; 

116 

x_max  =  xmax  >  x_max  ?  xmax  :  x_max; 

117 

y_min  =  ymin  <  y_min  ?  ymin  :  y_min; 

118 

y_max  =  ymax  >  y_max  ?  ymax  :  y_max; 

119 

> 

120 

const  double  DELTA_X  =  (  x_max  -  x_min  )  /  double(n_dim) ; 

121 

const  double  DELTA_Y  =  (  y_max  -  y_min  )  /  double(n_dim) ; 

122 

const  double  DELTA_A  =  DELTA_X  *  DELTA_Y ; 

123 

124 

va: : Vector  p,  py; 

125 

double  ap  =  0. ; 

126 

for  (  double  y  =  y_min;  y  <=  y_max;  y  +=  DELTA-X  )  { 
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127 

128  py  =  y  *  j ; 

129  for  (  double  x  =  x_min;  x  <=  x_max;  x  +=  DELTA_Y  )  { 

130 

131  p  =  x  *  i  +  py; 

132  for  (  register  int  n  =  0;  n  <  N_ FACETS;  n++  )  { 

133  vl  =  facet_vl[n]; 

134  v2  =  facet_v2[n]; 

135  v3  =  facet_v3[n]; 

136  vl=(vl*i)*i+(vl*j)*j; 

137  v2  =  (  v2  *  i  )  *  i  +  (  v2  *  j  )  *  j ; 

138  v3=(v3*i)*i+(v3*j  )*j; 

139  if  (  inside (  p,  vl,  v2,  m3  )  )  { 

140  ap  +=  DELTA-A; 

141  break; 

142  } 

143  } 

144  } 

145  } 

146  return  ap; 

147  } 

148 

149  int  main(  int  argc,  char*  argv[]  )  { 

150 

151  char  *  file; 

152  int  n_dim  =32;  //  grid  size  is  n_dim  x  n_dim,  with  32  x  32  =  1024  as  default 

153  if  (  argc  ==  3  )  { 

154  file  =  argvfl] ; 

155  n_dim  =  atoi(  argv[2]  ); 

156  > 

157  else  if  (  argc  ==  2  ) 

158  file  =  argv[l]; 

159  else  { 

160  std::cerr  «  "Usage:  "  «  argv[0]  «  "  stl_binary_inputfile"  «  std::endl; 

161  exit(  1  ); 

162  } 

163  char  header [100]  = 

164  unsigned  long  int  number  =  0; 

165  facet  tri; 

166 

167  std : : if stream  fin; 

168  fin. open (  file,  std::ios::in  |  std :: ios :: binary  ); 

169  if  (  fin. bad ()  )  { 

170  std::cerr  «  "Error  in  opening  file  "  «  file  «  std::endl; 

171  exit (  1  ); 

172  > 

173  std::cerr  «  "Contents  of  "  «  file  «  "  «  std::endl; 

174  fin.read(  (char  *)  header,  80  ); 

175  fin. read (  (char  *)  &number,  4  ); 

176  std:: clog  «  "Header  =  "  «  header  «  std::endl; 

177  std::clog  «  "There  are  "  «  number  «  "  triangles  in  "  «  file  «  ":  "  «  std::endl; 

178 

179  va:: Vector  n,  vl,  m2,  v3; 

180  double  surface_area  =  0.,  vol6  =  0.; 

181  std: :vector<va: :Vector>  facet_vl(number) ,  facet_v2 ( number ) ,  facet_v3 ( number) ; 

182 

183  for  (  unsigned  int  i  =  0;  i  <  number;  i++  )  { 

184 

185  fin. read (  (char  *)  &tri,  50  ); 

186 


187 

n  =  va : : Vector( 

(double)tri.nx, 

(double)tri.ny,  (double)tri.nz  ); 

188 

vl  =  va : : Vector( 

(double)tri . vlx, 

(double)tri.vly, 

(double)tri . viz  ); 

189 

m2  =  va : : Vector ( 

(double)tri.v2x, 

( double )tri.v2y, 

(double)tri. v2z  ); 

190 

v3  =  va : : Vector ( 

(double)tri.v3x, 

( double )tri.v3y, 

(double)tri.v3z  ); 

191 

192  surface_area  +=  area2(  n,  vl,  m2,  m3  ); 

193  vol6  +=  volume6(  vl,  m2,  m3  ) ; 

194 

195  facet_vl[i]  =  vl; 

1%  facet_v2[i]  =  v2; 

197  facet_v3[i]  =  v3; 

198  > 

199  f in. close( ); 

200  surface-area  /=  2.; 

201 

202  const  double  VOLUME  =  fabs(  vol6  /  6.  ); 

203  const  double  F  =  pow(  VOLUME,  -2.  /  3.  ); 

204 

205  std::cerr  «  "Volume  =  "  «  VOLUME  «  std::endl; 

206  std::cerr  «  "Number  of  facets  =  "  «  number  «  std::endl; 

207  std::cerr  «  "Total  surface  area  =  "  «  surface_area  «  std::endl; 

208  std::cerr  «  "Mean  shape  factor  (using  Cauchy's  theorem  for  convex  solid)  =  "  «  0.25  *  surface_area  *  F  «  std::endl; 

209 

210  const  int  N  =  10; 

211  double  x,  y,  z,  sf; 
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212 

213 

214 

215 

216 

217 

218 

219 

220 

221 


rng : : Random  rng; 

for  (  int  i  =  0;  i  <  N;  i++  )  { 


rng.spherical_avoidance(  x,  y,  z  ); 

sf  =  F  *  area_projected (  facet_vl,  facet_v2,  facet_v3,  va::Vector(  x,  y,  z  ),  n_dim  ); 
std::cout  «  sf  «  std::endl; 


return  EXIT-SUCCESS; 


Thus,  to  characterize  the  shape  factor  distribution  of  natural  fragments,  we  can 
either  measure  many  fragments  using  the  Icosahedron  Gage,  or  we  can  laser  scan, 
possibly  a  fewer  number.  Two  parameters  are  sufficient  to  characterize  the  lognormal 
distribution,  in  either  case.  Then  it  is  a  simple  matter  to  simulate  the  shape  factor 
using  Listing  6. 

4.  Shape  Factor  Modehng 

Finally,  let  us  consider  modeling  the  shape  factor.  Some  penetration  models,  such  as 
THOR,3  only  require  the  presented  area  of  the  fragment  at  impact,  in  which  case  the 
shape  factor  is  sufficient.  But  other  penetration  models,  such  as  FATEPEN,1’2  require 
a  specific  shape  and  orientation.  This  requires  a  realization  of  the  shape  factor,  and 
the  simplest  shapes  in  FATEPEN  that  allow  for  this  are  a  cylinder  and  a  cuboid. 

The  procedure  we  use  applies  to  both  shapes.  We  start  with  the  cylinder  and  cuboid 
in  standard  orientation,  which  is  with  the  axis  of  symmetry  along  the  2-axis  in  the 
case  of  the  cylinder,  and  with  the  length  along  the  2-axis,  width  along  the  x-axis, 
and  thickness  along  the  y- axis  in  the  case  of  the  cuboid.  The  target  is  taken  to  be  in 
the  x-y  plane.  Then  we  are  given  the  mass  m,  material  density  p,  and  shape  factor  7. 

For  cylinders,  we  randomly  select  a  candidate  Lj D,  which  then  allows  us  to  compute 
a  minimum  and  a  maximum  shape  factor  for  this  cylinder,  depending  upon  its 
orientation.  If  the  drawn  shape  factor  falls  within  this  range,  7min  <  7  <  7max,  then 
we  know  there  is  some  orientation  that  will  work;  in  the  next  section  we  work  out 
the  formula  for  the  appropriate  yaw  angle.  We  then  have  enough  information  to  also 
compute  the  diameter  and  the  length  of  the  cylinder. 

For  cuboids,  it  is  the  same  idea  except  that  we  need  to  select  candidate  W/L  and 
T jW  ratios,  which  then  provide  enough  information  to  compute  a  minimum  and 
a  maximum  shape  factor,  depending  upon  orientation.  We  then  compute  the  3D 
rotation  that  will  realize  this  shape  factor.  Finally,  we  factor  this  rotation  into  a 
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pitch-yaw-roll  rotation  sequence,13  as  these  are  required  for  FATEPEN. 

4.1  Cylinder 

The  formula  for  the  dimensionless  shape  factor  of  a  right-circular  cylinder  (RCC)  as 
a  function  of  the  effective  yaw  angle*  <j)y  is14 

=  asm  (j)y  +  b  cos  (f)y,  (19) 


where 


7 T  L 
OEI4  D 


-2/3 


D  ^  &E  4P 


7T  L\  2//'3  7T 


(20) 


The  yaw  angle  that  gives  the  maximum  shape  factor  is  obtained  by  setting  the 
derivative  with  respect  to  <py  equal  to  zero  and  solving  for  cj)y\ 


d'y 

/  -/= 


=  a  cos  4>y  —  b  sin  tj)y  =  0, 


which  gives 


7=7max 


'  V  7=7max 


=  tan 


-l 


To  simplify  the  notation,  let  oy  denote  this  angle:  (j)y  =  0j,7=7max.  Then, 

7max  =  7  {fy)  =  «  ,  ”,  ^  =  Vfl2  +  62, 

V«  +  b2  va  +  b2 

and  we  can  write 


(21) 


(22) 


(23) 


a  =  7max  sin  <j>y  and  b  =  ymax  cos  </>„, 


(24) 


so  that  Eq.  19  can  be  written  as 


Tmax  COs(0y  4>y)  if  4*y  $y 


'Tmax  COs((f)y  4*y)  if  4*y  ^ 


(25) 


*Note  that  effective  yaw  includes  pitch  and  is  defined  by  <f)yfin  =  cos  1(cos  <f>p  cos  (f>y),  where 
<f>p  is  pitch  and  <j>y  is  actual  yaw,  but  we  use  <f>y  rather  than  0y,eff  just  to  simplify  the  notation. 
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where  oy  =  cos  l(b/ 7max).  Solving  for  the  yaw  angle,  we  get 


(cos  1{b/% max) 

+  COS  1  (T  /  Tmax) 

if  7  <  6 

(t>y=\ 

^  COS  1(b/'f max) 

COS-1  (T  /  Tmax) 

IV 

cr* 

(26) 


which  shows  that  we  can  easily  get  the  orientation  (yaw  angle)  of  an  RCC  from  the 
shape  factor  at  impact. 

The  minimum  shape  factor  is 


Tmin  \ 


7T  L  5 

-2/s  L 

if 

L 

4  D  ) 

1  D 

Id 

7T  L  5 

-2/3 
|  7T 

if 

l 

4D ) 

1  4 

D 

7 r 
4 

7 r 

4 


(27) 


Solving  this  for  L/ D  for  a  given  minimum  shape  factor,  we  get 


L 

D 


L  7 r 

if  —  c  —  (disk-like) 


7T  \  X/2  _q/9  L  7 r 

)  7min  if  ^  >  4  (rod-like) 


(28) 


The  minimum  shape  factor  for  natural  fragments  is  around  0.5.  Inserting  this  value 
into  this  equation  gives 


L  7T 

r  ,  0.077  if  —  <  —  (disk-like) 

L  |  D  4 

D  |  L  7T 

2.5  if  —  >  —  (rod-like) 
4 


(29) 


The  maximum  shape  factor  for  natural  fragments  is  around  5.5.  Inserting  this  value 
into  the  appropriate  cubic  equation  (see  Appendix  A,  Eq.  A-46)  gives 


L  7 r 

,  0.0691055  if  <  —  (disk-like) 

Ju  |  LJ  4 

D  \  L  7T 

102.619  if  —  >  -  (rod-like) 
D  4 


(30) 


The  L/ D  value  for  a  rod- like  RCC  is  unrealistic,  so  we  choose  disk-like  RCCs  to 
model  the  maximum  shape  factor.  Selecting  disk-like  RCCs  with  L/D  in  the  range 
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from  0.069  to  7t/4  will  span  the  range  of  shape  factors  from  0.5  to  5.5.  The  actual 
code  that  generates  a  yawed  RCC  to  represent  the  shape  factor  is  given  in  Listing  12. 


Listing  12.  sf-rcc.cpp 

//  sf-rcc.cpp:  Implementation  of  an  algorithm  for  generating  FATEPEN  RCCs  to  represent  a  specified  shape  factor. 

//  Given  a  dimensionless  shape  factor,  generates  the  L/D  and  yaw  angle  for  the  RCC  to  represent  it. 

//  The  RCCs  are  disk-like  in  order  to  span  the  range  of  shape  factors  from  0.5  to  5.5. 

//  R.  Saucier,  October  2011 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

int  main(  int  argc,  char*  argv[]  )  { 

int  N  =  1000;  //  number  of  samples  or  override  on  command  line 

const  double  R2D  =  180.  /  M_PI;  //  to  convert  from  radians  to  degrees 

const  double  SF_MIN  =  0.5;  //  minimum  shape  factor  (found  from  artillery  fragments) 

const  double  SF_MAX  =  4.5;  //  maximum  shape  factor  (found  from  artillery  fragments) 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  );  //  override  default  number  of  samples  on  command  line 

//  default  values  for  the  shape  factor  lognormal  distribution  from  122mm,  152mm  and  155mm  artillery 
double  mu  =  0.596514;  //  these  two  parameters  characterize  the  lognormal  shape  factor  distribution 

double  sigma  =  0.340874;  //  with  mode  =  1.62,  median  =  1.81  and  mean  =  1.93 

rng : : Random  rng; 

double  a,  b,  c,  th,  sf_min,  sf_max,  sf,  l_d,  l,  d,  yaw,  V  =  1.;  //  here  we  use  a  fragment  with  unit  volume 

for  (  int  n  =  0;  n  <  N;  n++  )  { 

//  normally,  the  shape  factor  would  be  provided,  but  here  we  get  a  shape  factor  within  bounds  [SF_MIN,  SF_MAX] 
do  {  sf  =  rng. lognormal (  0.,  mu,  sigma  );  }  while  (  sf  <  SF_MIN  | |  sf  >  SF-MAX  ); 

//  now  we  want  to  realize  this  shape  factor  with  a  yawed  cylinder 
do  { 

l_d  =  rng.uniform(  0.069,  M_PI_4  );  //  disk-like  RCCs  (l_d  =  0.069  corresponds  to  sf_max  =  5.5) 
c  =  pow(  M_PI_4  *  l_d,  -2./3.  ); 
a  =  c  *  l_d; 
b  =  c  *  M_PI_4; 
sf_min  =  a; 

sf_max  =sqrt(  a*a+b*b  ); 

}  while  (  sf  <  sf_min  | |  sf  >  sf_max  ) ; 

if  (  sf  <  b  )  th  =  acos(  b  /  sf_max  )  +  acos(  sf  /  sf_max  ); 
else  th  =  acos(  b  /  sf_max  )  -  acos(  sf  /  sf_max  ); 

d  =  pow (  V  /  (  M_PI_4  *  l_d  ),  1./3.  ); 
l  =  d  *  l_d; 
yaw  =  th  *  R2D; 

std::cout  «  sf  «  "\t"  «  l_d  «  "\t"  «  yaw  «  std::endl; 

} 

return  EXIT.SUCCESS; 


This  algorithm  is  stochastic,  so  that  for  a  given  shape  factor  there  will  be  a  range 
of  L/D  ratios  and  yaw  angles  that  can  represent  the  fragment.  This  will  result  in  a 
range  of  residual  velocities,  but  this  range  will  be  relatively  small.  For  example,  a 
725-gr  steel  RCC  with  a  shape  factor  of  1.93  striking  a  1/4-inch  mild  steel  plate  at 
3500  f/s  results  in  a  residual  velocity  of  1988  ±  23  f/s. 

4.2  Cuboid 

The  procedure  for  realizing  the  shape  factor  as  a  cuboid  with  a  specific  pitch,  yaw, 
and  roll15  is  implemented  in  the  code  in  Listing  13. 
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Listing  13.  sf-rpp.cpp 


1  //  sf-rpp.cpp:  Monte  Carlo  shape  factor  from  a  lognormal  distribution 

2  //  This  code  demonstrates  how  it  is  possible  to  make  use  of  a  lognormal  shape  factor  distribution, 

3  //  and  at  the  same  time  create  RPPs  for  FATEPEN  that  have  the  proper  volume  and  presented  area. 

4  //  Given  the  fragment  volume  (or  mass  and  density),  it  provides  everything  that  FATEPEN  requires: 

5  //  the  length,  width,  and  thickness  of  the  RPP,  as  well  as  the  pitch -yaw- roll  rotation  sequence 

6  //  that  will  take  the  RPP  from  standard  orientation  to  the  orientation  that  realizes  the  shape  factor. 

7  //  R.  Saucier,  September  2011 

8 

9  #include  "Random. h" 

10  #include  "Rotation . h" 

11  #include  <iostream> 

12  #include  <cstdlib> 

13  #include  <cassert> 

14  #include  <iomanip> 

15 

16  const  double  SF_MIN  =  0.5;  //  minimum  shape  factor 

17  const  double  SF_MAX  =  4.5;  //  maximum  shape  factor 

18  const  double  MU  =  0.596514;  //  these  two  parameters  characterize  the  lognormal  shape  factor  distribution 

19  const  double  SIGMA  =  0.340874;  //  mode  =  1.62,  median  =  1.81,  mean  =  1.93 

20  const  va : : Vector  I(  1.,  0.,  0.  ),  J(  0.,  1.,  0.  ),  K(  0.,  0.,  1.  ); 

21 

22  int  main(  void  )  { 

23 

24  const  int  N_SAMPLES 

25  const  double  W_L_MIN 

26  const  double  T_W_MIN 

27  const  double  G2GR 

28  const  double  GR2G 

29  const  double  MASS 

30  const  double  RHO 

31  const  double  V 

32  const  double  V_23 

33 

34  double  th,  th_max,  Amin,  Amax,  w_l,  t_w.  Ax,  Ay,  Az,  W,  L,  T; 

35  va:: Vector  Ap,  u_max,  axis,  u; 

36  va:: Rotation  R; 

37  va:: sequence  s; 

38  rng:: Random  rng; 

39  double  cp,  cy,  cr,  sp,  sy,  sr,  p,  y,  r,  sf,  ap,  apcalc; 

40 

41  std::cout  «  std : : setprecision (6)  «  std:: fixed; 

42  for  (  int  n  =  0;  n  <  N_SAMPLES;  n++  )  { 

43 

44  //  normally,  the  shape  factor  would  be  provided,  but  here  we  get  a  shape  factor  within  bounds  [SF-MIN,  SF-MAX] 

45  do  {  sf  =  rng .lognormal (  0.,  MU,  SIGMA  );  }  while  (  sf  <  SF_MIN  | |  sf  >  SF_MAX  ); 

46 

47  sf  =  1.93;  //  or  select  from  a  lognormal  distribution 

48  ap  =  sf  *  V_23; 

49 

50  do  { 

51  w_l  =  rng.uniform(  W_L_MIN,  1.  ); 

52  t_w  =  rng.uniform(  T_W_MIN,  1.  ); 

53  W  pow(  V  *  w_l  /  t_w,  1./3.  ); 

54  L  =  W  /  w_l; 

55  T  =  W  *  t_w; 

56  Ax  =  L  *  T; 

orientation 

57  Ay  =  L  *  W; 

58  Az  =  W  *  T; 

59  Amin  =  Az; 

60  Amax  =  sqrt(  Ax  *  Ax  +  Ay  *  Ay  +  Az  *  Az  ) 

61  }  while  (  ap  <  Amin  | |  ap  >  Amax  ) ; 

62 

63  Ap  =  Ax  *  I  +  Ay  *  J  +  Az  *  K;  //  presented  area  vector 

64  U-max  =  Ap  /  Amax;  //  u_max  is  direction  which  realizes  Amax  =  Ap  *  u_max; 

65  R  va :: Rotation (  K,  U-max  );  //  rotation  which  takes  K  to  u_max,  and  will  take  |Ap|  from  Amin  to 

Amax 

66  axis  =  va: : Vector (  R  );  //  fixed  axis  of  this  rotation 

67  th_max  =  double (  R  );  //  angle  of  rotation  to  rotate  K  to  u_max 

68  th  th_max  -  acos(  ap  /  Amax  );  //  rotation  angle  for  u 

69  R  va: : Rotation (  axis,  th  ) ;  //  rotation  which  takes  K  to  u 

70  va::Rotation  R2(  K,  rng . uniform(  0.,  2.  *  M_PI  )  );  //  this  should  have  no  effect  upon  the  projected  area  perp  to 

K 


//  ratio  of  W./L 
//  ratio  of  T/W 
//  width  of  RPP 
//  length  of  RPP 
//  thickness  of  RPP 

//  presented  area  orthogonal  to  x-axis  (intermediate)  in  initial 

//  presented  area  orthogonal  to  y-axis  (maximum) 

//  presented  area  orthogonal  to  z-axis  (minimum) 

//  min  presented  area 
//  max  presented  area 


=  10000; 

=  0.185; 

// 

minimum  W/L 

ratio 

(must 

be  0.185 

or  smaller 

to 

get 

SF_MAX  = 

5.5) 

=  0.185; 

// 

minimum  T/W 

ratio 

(must 

be  0.185 

or  smaller 

to 

get 

SF_MAX  = 

5.5) 

=  15.4324;  //  to  convert  grams  to  grains 

=  1.  /  G2GR;  //  to  convert  grains  to  grams 

=  725.  *  GR2G;  //  mass  (725  grains  converted  to  grams) 

=  7.83;  //  density  of  steel  (g/cirT3) 

=  MASS  /  RHO;  //  volume  (cmA3) 

=  pow (  V,  2./3.  ); 


71  s  va::factor(  R2  *  R,  va::XYZ  );  //  pitch -yaw- roll  rotation  sequence 

72 

73  p  =  s. first;  //  pitch  (rad) 

74  y  =  s. second;  //  yaw  (rad) 

75  r  =  s. third;  //  roll  (rad) 


76 

cp  =  cos( 

p  ); 

77 

cy  =  cos( 

y  ); 

78 

cr  =  cos( 

r  ); 

79 

sp  =  sin( 

p  ); 

80 

sy  =  sin( 

y  ); 
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81 

sr  =  sin(  r  ) ; 

82 

apcalc  =  fabs(  sp  * 

sr  -  cp 

*  sy  *  cr  ) 

*  Ax  +  fabs(  sp  * 

cr  +  cp  *  sy  *  sr  )  *  Ay  +  fabs(  cp  *  cy  )  *  Az; 

83 

assert  (  fabs(  ap  - 

apcalc  ; 

1  <  0.001  ); 

84 

std : : cout  «  w_l  « 

"\t"  « 

t_w  «  "\t" 

«  p  *  va : :  R2D  « 

"\t"  «  y  *  va::R2D  «  "\t"  «  r  *  va::R2D  «  std: 

85 

86 

> 

87 

return  EXIT-SUCCESS; 

88 

} 

This  will  result  in  a  range  of  residual  velocities.  For  example,  a  725-gr  steel  cuboid 
with  a  shape  factor  of  1.93  striking  a  1/4-inch  mild  steel  plate  at  3500  f/s  results 
in  a  residual  velocity  of  2464  ±  363  f/s,  much  more  variation  than  is  the  case  with 
cylinders. 

5.  Conclusions  and  Recommendations 

We  have  provided  explicit  analytical  formulas  for  the  shape  factor  distributions  of 
some  common  shapes  with  random  orientations.  And  we  have  shown  that  it  is  easy 
to  simulate  these  shape  factor  distributions  with  computer  code  and  demonstrated 
through  plots  that  the  simulations  match  the  plots  from  the  analytical  formulas.  When 
we  examine  natural  fragment  shape  factor  distributions,  however,  the  only  shape  that 
comes  close  is  an  ellipsoid.  But  it,  too,  fails  to  provide  an  adequate  representation  by 
simply  randomizing  its  orientation.  We  can  work  out  the  dimensions  of  the  ellipsoid 
from  the  projected  area  measurements,  but  then  the  volume  comes  out  wrong  because 
it  does  not  account  for  hidden  surfaces. 

A  better  approach  to  shape  factor  simulation  was  found  after  enough  natural  fragment 
data  was  processed  to  reveal  that  it  could  be  fit  with  a  lognormal  distribution.  The 
measurement  of  fragment  shape  factor  with  the  Icosahedron  Gage  does  not  give 
any  indication  of  a  lognormal  distribution  with  only  16  measurements,  but  once 
we  combine  the  measurements  from  hundreds  of  fragments,  the  resemblance  to 
the  lognormal  is  striking.  Rather  than  trying  to  find  a  shape  that  will  will  work 
by  randomizing  the  orientation,  it  makes  more  sense  to  use  the  lognormal  as  a 
probability  distribution  in  Monte  Carlo  sampling. 

We  also  showed  that  laser  scans  of  fragments  can  be  used  to  compute  the  fragment 
shape  factor  from  any  viewpoint,  and  we  described  a  variety  of  methods  of  achieving 
a  uniform  spherical  distribution.  Computing  fragment  volume  and  projected  area  is 
fast  if  we  treat  the  fragments  as  convex  solids.  But  natural  fragments  are  not  convex, 
and  accounting  for  hidden  surfaces  uses  much  more  computer  time. 
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Finally,  we  showed  that  it  is  possible  to  realize  each  fragment  mass  and  shape  factor 
as  either  a  yawed  cylinder  or  a  cuboid  with  a  pitch,  yaw,  and  roll.  Thus,  we  have  a 
procedure  for  generating  all  the  input  variables  required  to  run  THOR  or  FATEPEN 
from  the  given  mass  and  shape  factor. 
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Appendix  A.  Analytical  Shape  Factor  Formulas  for  5  Convex  Solids 
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Here  we  cite  explicit  formulas  for  the  probability  density  function  (PDF)  and  the 
cumulative  distribution  function  (CDF)  for  5  convex  solids  that  have  a  random 
orientation.  For  each  solid,  we  also  list  code  that  can  be  used  to  plot  the  PDF  and 
CDF 


A-1.  Cube 

The  PDF  /( 7)  is  given  by1 


/( 7)  = 


4^3 

- tan 

7 r 


'V3 


F(7)  =  < 


-  6 


tan 


7 r 

'2- vV 


il-3 

73 


7 


7^ 


and  the  CDF  F( 7)  is  given  by 


47  4737  _! 

— j=  —  3 - tan 

73 


^7^7 

7 


if  1  <  7  <  72 
if  72  <  7  <  73 


tan 


2  +  737 

.v/2^7 


if  1  <  7  <  72 
if  72  <  7  <  73 


(A-1) 


(A-2) 


These  formulas  are  implemented  in  the  code  in  Listing  A-1  and  shown  plotted  in 
Fig.  A-1. 


Listing  A-1.  cube.cpp 


//  cube.cpp:  generates  plotting  points  for  cube  pdf  and  cdf 
//  Ref:  Vickers,  G.  T.  and  Brown,  D.  J., 

//  "The  distribution  of  projected  area  and  perimeter  of  convex,  solid  particles," 

//  Proc.  R.  Soc.  Lond.  A  (2001),  Vol.  457,  pp.  283-306. 

#include  <iostream> 

#include  <cmath> 

#include  <cstdlib> 

#include  <cassert> 

static  const  double  SQRT2  =  M_SQRT2,  SQRT3  =  sqrt(  3.  ); 

double  pdf (  double  x  )  { 

assert (  1.  <=  x  &&  x  <=  SQRT3  ); 

if  (  x  <=  SQRT2  )  { 

double  a  =  sqrt(  2.  -  x  *  x  ); 

return  4.  /  SQRT3  -  4.  *  SQRT3  *  atan(  SQRT3  *  a  /  x  )  /  M_PI; 

> 

else 

return  4.  /  SQRT3 ; 

} 

double  cdf(  double  x  )  { 

assert (  1.  <=  x  &&  x  <=  SQRT3  ); 

if  (  x  <  SQRT2  )  { 

double  a  =  sqrt(  2.  -  x  *  x  ); 


'Vickers  GT,  Brown  DJ.  The  distribution  of  projected  area  and  perimeter  of  convex,  solid  particles. 
Proc  R  Soc  Lond  A.  2001;457:283-306. 
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31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 


return  4.  *  x  /  SQRT3  -  3.  -  (  4.  *  SQRT3  *  x  /  M_PI  )  *  atan(  SQRT3  *  a  /  x  )  + 

(  6.  /  M_PI  )  *  (  atan(  (  2.  -  SQRT3  *  x  )  /  a  )  +  atan(  (  2.  +  SQRT3  *  x  )  /  a  )  ); 

} 

else 

return  4.  *  x  /  SQRT3  -  3.; 


int  main(  int  argc,  char*  argv[]  )  {  //  specify  number  of  points  on  commandline  or  use  1000 

int  N  =  1000; 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  ); 
const  double  GMIN  =1.,  GMAX  =  SQRT3 ; 

for  (  double  g  =  GMIN;  g  <=  GMAX;  g  +=  (  GMAX  -  GMIN  )  /  N  ) 

std::cout  «  g  «  "\t"  «  pdf(  g  )  «  "\t"  «  cdf(  g  )  «  std::endl; 

return  EXIT_SUCCESS; 

> 


The  conventional  way  of  turning  this  into  a  shape  factor  probability  distribution  is 
to  sample  the  uniform  distribution  between  0  and  1,  P  ~  [7(0, 1),  and  then  invert 
F  to  get  7  =  F~]  (P).  But  there  is  no  simple  way  to  solve  Eq.  A-2  for  7  when 
1  <  7  <  v/2,  so  instead  we  simulate  the  shape  factor  probability  distribution  by 
uniform  random  sampling  over  the  unit  sphere  and  compute  the  shape  factor  for 
each  orientation  to  build  up  a  probability  distribution. 


Fig.  A-l.  Plot  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  cube 


A-2.  Cuboid 

The  source  for  the  following  formulas  is  Walters.2 

A-2.1.  Notation 

Let  L  >  W  >  T  be  the  length,  width,  and  thickness  of  the  cuboid,  so  that  length  is 
always  the  largest  dimension  and  thickness  the  smallest.  By  setting 

Ax  =  WT,  Ay  =  TL,  Az  =  LW,  (A-3) 

2Walters  AG.  The  distribution  of  projected  areas  of  fragments.  Proc  Cambridge  Phil  Soc  Math 
and  Phys  Sci.  1947;43:343-347. 
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we  order  the  areas  of  the  3  faces  so  that  Ax  <  Ay  <  Az  and  match  the  notation  in 
Walters2.  If  we  only  know  the  3  face  areas,  then  we  can  get  the  dimensions  from 


L 


AyAz 


Ar 


W 


AXAy 


A, 


(A-4) 


To  simplify  the  notation,  let  x  represent  the  variable  presented  area  and  define 


a  =  Ax,  b  =  Av ,  c  =  Az,  nr  =  a2  +  b2  +  c2, 


(A- 5) 


x. 


=  \/m2  —  a2,  Xb  =  \/m 2  —  62,  xc  =  '/m2  —  c2. 


(A- 6) 


Walters2  also  defines  the  following  quantities: 


2i'”2  ™2_l,2  2  m2x2 


_  m  -\-  x  a  —  m  +  c' 
^  =  .2  _  T2  ’  ^  = 


rrr  —  ar 


m2  —  c2  ’ 


- 2’  D*  = 

rnz  —  xz 


2  m2c2 
in'2  —  c2  ’ 


Q2!  = 


C2  +  X2 


«T  = 


9  i  95 

c2  +  a2 


c2  +  62 


(A-7) 

(A- 8) 


A-2.2.  Integrals 

Using  x  as  the  variable  presented  area,  we  make  use  of  the  following  4  indefinite 
integrals: 


h(x)  =  j  sin  1(A  —  aiD)dx 
=  a:  sin-1  (A  —  a.\ D) 


+  2 to  tan 


-1  (  'Ac  -  x2 


—  2ctanh 


-1 


\/  x2c  -  x2 

TO 


for  a  <  x  <  xc  (A-9) 


I2(x)  =  J  sin  1(A  —  axD)dx 

=  xsin_1(A  —  axD )  —  2?ntan 

I3(x)  =  J  sin_1(A  —  ayD)dx 

=  xsin~  (A  —  ayD)  —  2?ntan 

litx)  =  J  sin~1(Az  —  a\ Dz)dx 

=  xsin^1(Az  —  ot\Dz)  +  2ctanh 


-1  [  vA  - 2:2 


for  a  <  x  <  Xb 


-1  (  Vx2a~' 


for  a  <  x  <  xa 


vA  - 


for  a  <  x  <  xc 


(A- 10) 


(A- 11) 


TO- 


(A- 12) 
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The  corresponding  definite  integrals  are 


x2)  =  h{x 2)  -  h{x  1)  for  i  =  1,  2, 3, 4.  (A- 13) 

It  is  also  convenient  to  define  the  following  2  constants: 

k\  =  sin-1  (A.  —  axDz)  and  k2  =  sin_1(Az  —  ayDz).  (A- 14) 

A-2.3.  Probability  Density  Function  and  Cumulative  Distribution  Function 

Let  f(x)  represent  the  probability  density  function  and  Fix)  represent  the  cumulative 
distribution  function.  The  following  expressions  for  f{x)  are  taken  from  Walters,2 
and  F(x)  is  obtained  by  integrating  the  corresponding  density  function.  There  are  2 
cases  to  consider,  depending  upon  the  value  of  a 2  +  b 2  relative  to  c2,  and  each  case 
has  6  distinct  regions. 

Case  1:  a2  +  b2  <  c2 

•  Case  1.1:  a  <  x  <  b 

f(x)  =  — [sin_1(A  —  ot\ D)  —  sin_1(A  —  axD)  +  sin_1(A,  —  a\Dz)  —  k\\  (A-15) 

7TTO 

F{x)  =  - [Ii(a,x)  —  I2(a,x)  +  I4  {a,x)  —  k\{x  —  a)]  (A-16) 

7 rra 

•  Case  1.2:  b  <  x  <  xc 

}{x)  =  ^*—[2  sin-1  (A  —  ai  D)  —  sin_1(A  —  axD)— 
nm 

sin-1  (A  —  ayD)  +  2sin_1(Az  —  ot\ Dz)  —  k\  —  k2] 

F(x)  =  —  [h{a,  b )  -  J2(a,  b)  +  I^a,  b)  —  ki(b  —  a)  +  2I1(b,x)- 
nm 

I2(b,  x)  —  I3 (6,  x)  +  2/4(6,  x)  -  ki(x  —  b)  —  k2(x  —  6)] 

•  Case  1.3:  xc  <  x  <  c 

f(x)  =  — [2n  —  sin_1(A  —  axD)  —  sin_1(A  —  avD)  —  k\  —  k2]  (A-19) 

nm 

F{x)  = - [h(a,  b)  -  I2(a,  b)  +  h(a,  b)  —  k\{b  -  a)  +  2/i(6,  xc)— 

nm 

h{b,xc)  -  I3(b,  xc)  +  2/4(6,  xc)  -  k\(xc  -  6)  -  k2(xc  -  b)+ 

2n(x  -  xc)  -  I2(xc,  x)  -  h(xc,  x)  -  k\{x  -  xc)  -  k2{x  -  xc)\  (A-20) 


(A- 17) 

(A- 18) 
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•  Case  1.4:  c  <  x  <  xb 


2 

fix)  =  - [7 r  —  sin_1(A  —  axD)  —  sin_1(A  —  avD)\ 

7TTO 

F(x)  =  —  [h(a,  6)  -  I2 (a,b)  +  h(a,b)  —  ki(b  — a)  +  2I1(b,xc)— 
7 xm 

h{b,  xc)  -  h(b,  xc)  +  2/4(6,  xc)  -  k i(xc  -  6)  -  fc2(xc 
27 r(c  -  xc)  -  h(xc,  c)  -  h{xc,  c)  -  fci(c  -  xc)  -  k2(c 
2n(x  —  c)  —  2/2(0,  x)  —  2/3(0,  x)] 

•  Case  1.5:  xb  <  x  <  xa 

fix)  =  —  [37t  —  2  sin-1  (A  —  avD)] 
nm 

F{x)  = - [/i(a,  6)  —  /2(a,  6)  +  h(a,b)  —  £4(6  —  a)  +  2/i(6,xc)— 

7r  m 

hip,  xc)  -  h(b,  xc)  +  2/4(6,  xc)  -  ki(xc  -  6)  -  k2{xc 
27 r(c  -  xc)  -  h{xc,  c)  -  /3(®c,  c)  -  fci(c  -  xc)  -  fc2(c 
27r(xft  -  c)  -  2/2 (c,  xfc)  -  2/3(c,  xb)+ 

37t(x  -  xb)  -  2/3 (xb,  x)] 


•  Case  1.6:  xa  <  x  <  m 


F{x)  =  —  [/i(a,  6)  -  /2(a,  6)  +  h{a,b)  -  £4(6  -  a)  +  2/i(6,  xc)— 
nm 

hip,  xc)  -  h(b,  xc)  +  2/4(6,  xc)  -  fci(xc  -  6)  -  fc2(xc 
27t(c  -  xc)  -  I2{xc,  c)  -  /3(®c,  c)  -  £4(0  -  xc)  -  fc2(c 
27r(xb  -  c)  -  2/2(c,  Xb)  -  2/3(c,  xb)+ 

4 

37t(xq  -  Xb)  -  2/3(Xb,Xa)]  H - (x  -  Xa) 


Case  2:  a2  +  b2  >  c2 

•  Case  2.1:  a  <  x  <  b 

fix)  = - [sin_1(A  —  ai-D)  —  sin_1(A  —  axD)  +  sin_1(A2  —  ot\Dz) 

nm 

F{x)  = - [/i(a,  x)  —  /2(a,  x)  +  h(a,x)  —  k\(x  —  a)] 

nm 
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b)+ 

xc)+ 

(A-22) 


(A-23) 

&)+ 

xc)+ 

(A-24) 


(A-25) 

b)+ 

xc)+ 

(A-26) 


£4]  (A-27) 
(A-28) 
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•  Case  2.2:  b  <  x  <  c 


/(x)  =  ——[2 sin  X(A  —  a\D)  —  sin  1  (A  —  axD)  —  sin  VA  —  avD)+ 
nm 

2sin-1(Az  —  ai  Dz)  —  k\  —  k2\ 

F(x)  =  —  [h{a,  6)  -  J2(a,  b )  +  /4(a,  6)  —  ki(b  —  a)  +  2/4(6,  x)— 

77771 

I2{b,  x)  -  I3(b,  x)  +  2/4(6,  x)  -  k\{x  -  b)  -  k2(x  -  6)] 

•  Case  2.3:  c  <  x  <  xc 

2 

f(x)  =  - [sin-1  (A  —  ol\ D)  —  sin-1  (A  —  axD)+ 

7 rm 

sin-1(Az  —  «i/)z)  —  sin-1  (A  —  ctyD)] 

F(x)  =  - [/i(a,  6)  —  /2(a,  6)  +  I^a,  b)  —  k\(b  —  a)  +  2/i(6,  c)— 

7 rm 

/2(6,  c)  -  /3(6,  c)  +  2/4(6,  c)  —  fci(c  —  6)  —  k2(c  -  b)+ 
2/4(0,  x)  -  2 12(c,x)  -  2/3 (c,  x)  +  2/4 (c,  x)] 

•  Case  2.4:  xc  <  x  <  xb 

2 

/(x)  =  - [77  —  sin-1  (A  —  a*/?)  —  sin-1  (A  —  a„D)l 

77777 

F(x)  =  [/i(a,6)  -  /2(a,  6)  +  /4(a,  6)  -  fci(6  -  a)  +  2/4(6,  c)— 

77777 

I2(b,  c)  -  /3 (6,  c)  +  2/4(6,  c)  -  k± (c  6)  -  k2{c~  6)+ 
2/4(0,  xc)  -  2/2(c,  xc)  -  2/3(c,  xc)  +  2/4(c,  xc)+ 

27r(x  -  xc)  -  2 /2(xc,  x)  -  2 /3(xc, x)] 

•  Case  2.5:  xb  <  x  <  xa 

f(x)  =  — - —  [877  —  2  sin-1  (A  —  a„D)l 

77  TO 

F(x)  = - [/i(a,  6)  —  /2(a,  6)  +  I^a,  b)  —  k\(b  —  a)  +  2/4(6,  c)— 

77  777 

/2(6,  c)  -  /3 (6,  c)  +  2/4(6,  c)  -  ki(c-b)  -  k2(c  -  b)+ 
2/4(0, xc)  -  2/2(c, xc)  -  2/3(o,xc)  +  2/4(c, xc)+ 

277(xb  -  xc)  -  2 /2(xc,  x6)  -  2/3 (xc,  x6)  + 

377(x  -  X6)  -  2/3 (Xfc,  x)] 
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•  Case  2.6:  xa  <  x  <  m 


f(x)  =  -  (A-37) 

m 

F(x )  =  - [/i(o,  6)  -  J2(a,  6)  +  14(0,  6)  —  k\(b  -  a)  +  2/i(6,  c)  — 

nm 

I2(b ,  c)  -  J3 (6,  c)  +  2/4(6,  c)  —  fci(c  —  6)  —  k2(c  -  b)+ 

2/i (c,  xc)  -  2/2(c,  xc)  -  2 13(c,xc)  +  2/4 (c,  xc)+ 

2ir(xb  -  xc)  -  2 I2(xc,  Xb)  -  2/3 (xc,  Xb)+ 

4 

37r(xa  -  xft)  -  2/3(xh,xa)]  4 - (x  -  xa)  (A-38) 

171 

These  formulas  have  been  implemented  in  Listing  A-2.  Example  plots  are  displayed 
in  Fig.  A-2. 


Fig.  A-2.  Plot  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  cuboid  with  L  —  3, 

FL  =  2,  T  =  1 
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//  rpp.cpp:  Given  the  side  lengths  of  an  RPP,  computes  the  shape  factor 
//  cumulativce  distribution  function,  which  is  is  obtained  by 

//  integrating  Walters'  formula  for  the  probability  density  function. 
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inline  double  i2(  double  a,  double  b,  double  c,  double  xl,  double  x2  )  {  return  i2( 


double  i3(  double  a,  double  b,  double  c,  double 
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:  3)  8  S  8 


double  a  =  W  *  T;  //  smallest  area,  Ax  in  Walters'  formula 

double  b  =  T  *  L;  //  intermediate  area,  Ay  in  Walters'  formula 


double  c  =  L  *  W;  //  largest  area,  Az  in  Walters'  formula 
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else  if  (  b  <=  x  &&  x  <  xc  )  {  //  case 


else  if  (  xc  <=  x  &&  x  <  c  )  {  //  case  c 

pdf  =  C  *  (  2.  *  M_PI  -  f2  -  f3  -  K1  -  K2  ); 
cdf  =C1+C2+C*  (  2.  *  M_PI  *  (  x  -  xc  )  -  i2(  a,  b,  c,  xc, 


x 


X 


(D 

'd' 


X 


ro 

ro 


m  o  h  oo  ! 


3  \o  ' 


lO'O^MW'tin'OC'OOO'i 


|  oo  oo  oo  ( 


;  VO  I 

i  In  OV  OV  I 


is  §  §  i 
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sf  =  x  *  S2;  //  convert  presented  area  to  dimensionless  shape  factor 

std::cout  «  sf  «  n\t"  «  pdf  *  SI  «  "\t"  «  cdf  «  std::endl;  //  shape  factor  PDF  and  CDF 


A-3.  Cylinder 


Consider  a  right-circular  cylinder  (RCC)  with  length  L  and  diameter  D.  Let  oy  be 
the  yaw  angle  measured  from  the  axis  of  symmetry  so  that  4>y  =  0  corresponds  to  a 
face-forward  orientation  of  the  cylinder.  Then  the  shape  factor  as  a  function  of  yaw 
angle  is3 


7  (02/) 


(A-39) 


where  0  <  <j)y  <  rr.  The  minimum  shape  factor  is 


7min 


-2/3 

min 


and  is  realized  at  the  yaw  angle 


<\>y 


7t/2  if  L/D  <  7t/4 
0  if  L/D  >  7t/4 


(A-40) 


(A-41) 


The  maximum  shape  factor  is 


(A-42) 


and  is  realized  at 

,  -i  (L/D\ 

The  mean  shape  factor  when  averaged  over  all  random  orientations  is 


(A-43) 


(A-44) 


Some  plots  of  Eqs.  A-39,  A-42,  and  A-44  as  a  function  of  L/D  are  shown  in 
Fig.  A-3. 


3Saucier  R.  Shape  factor  of  a  randomly  oriented  cylinder.  Aberdeen  Proving  Ground  (MD):  Army 
Research  Laboratory  (US);  2000  Jul.  Report  No.:  ARL-TR-2269. 
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L/D 


—  Maximum 

—  60  deg 

—  Mean 

—  45  deg 

—  30  deg 


Fig.  A-3.  Shape  factors  of  a  cylinder  as  a  function  of  L/D  at  a  fixed  yaw  angle  from  Eq.  A-39 
are  displayed  at  30°,  45°,  and  60°.  The  maximum  shape  factor  as  a  function  of  L/D  is  from 
Eq.  A-42,  and  the  average  shape  factor  of  a  randomly  oriented  cylinder  as  a  function  of  L/D 
is  from  Eq.  A-44. 


-I  /  Q 

The  minimum  of  the  maximum  shape  factor  curve  is  (7max)mm  =  y/3  (|)  ~ 

1.49  and  occurs  when  L/D  =  \/2  (f )  ~  1.11072.  The  minimum  of  the  mean  shape 
factor  curve  is  7min  =  |  (f )  ~  1.38395  and  occurs  when  L/D  =  1. 


In  Eq.  A-39,  let  x  =  (L/D)  1//3;  then  it  can  be  written  as 

"7T\-1/3  7 


X 


,  x  +  tan0„  =  0. 

cos  (py  V  4  / 


Square  Eq.  A-42  and  let  x  =  (L/D)  2/3;  then  it  can  be  written  as 

"K  \  — 2/3 


X  ~  .4 


7 r 


7max^'  +  (  J )  =  0. 


-2 


And  in  Eq.  A-44,  let  x  =  (L/D)  3/3;  then  it  can  be  written  as 


x6  -  2  - 


71 


-1/3 


'yx  +  2  =  0. 


So  we  see  that  these  3  equations  all  have  the  same  form: 


(A-45) 


(A-46) 


(A-47) 


x3  +  px  +  q  =  0, 


(A-48) 
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where 


x  —  — 


L 


-1/3 


x=|- 


D  / 

L  \  -2/3 


,P  =  - 


7T 


-1/3 


-1 


47  tan^inEq'A'45; 


>P  =  ~ 


7T  \  — 2/3 


7T 


-2 


Tmax’  ?  =  T  in  Eq-  A'46; 


X  =  — 


-1/3 


/7T\-V3  _ 

,p  —  —  2  J  7,  q  =  2  in  Eq.  A-47. 


A-3.1.  Diversion:  Solution  to  Cubic  Equation  in  the  Case  of  Real  Roots 

Without  loss  of  generality,  the  general  cubic  can  be  written  as 

a:3  +  ax 2  +  bx  +  c  =  0.  (A-49) 

quadratic  term  and  puts  it  in  the  form 
+  py  +  q  =  0,  (A-50) 

,  2 a3  .  _ . . 

and  g  =  c-y  +  — .  (A-51) 

We  note  in  passing  that  the  absence  of  the  quadratic  term  implies  that  the  roots  must 
sum  to  zero.  For  if  yu  y2,  and  y3  are  the  roots  of  Eq.  A-50,  then 


Setting  x  =  y  —  a/ 3  eliminates  the 


r 


with 


P  =  h~j 


{y-y\){y-y*){y-yz)  =  2/3 — (2/1 +2/2 +2/3)2/2 + (2/12/2 +2/12/3+2/22/3)2/ —2/12/22/3  =  0, 

(A-52) 


and  eliminating  the  y2  term  means  that  y1  +  y2  +  y3  —  0. 


A  useful  trick4  for  solving  the  cubic  without  the  quadratic  term  (Eq.  A-50)  is  to  make 
use  of  the  trigonometric  identity 


4  cos3  6  —  3  cos  9  —  cos(3 9)  =  0 


(A-53) 


which  can  be  easily  derived  by  using  the  double  angle  formulas  in  the  expansion  of 
cos(30).  Let  us  return  then  to  Eq.  A-48  to  see  if  we  can  transform  it  into  this  form. 

4Hubbard  JH,  Hubbard  BB.  Vector  calculus,  linear  algebra,  and  differential  forms:  a  unified 
approach.  Upper  Saddle  River  (NJ):  Prentice  Hall;  2001. 
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Begin  by  setting 


x  =  tcosO, 


(A-54) 


which  gives 

t3  cos3  9  +  pt  cos  9  +  q  =  0, 

or,  multiplying  through  by  4/t3, 


(A-55) 


4  cos3  9  +  -^  cos  9  +  ^  =  0. 
t 2  t6 


Choosing 


then  gives 


,  ,  4 p 

t  =  xl~Y 


4  cos3  9  —  3  cos  9  =  0. 


p  \  Ap 


(A-56) 


(A-57) 


(A-58) 


Now  we  see  that  if  we  choose  9  such  that 


COS(30)  =  :]!\-Ap 


(A-59) 


then  the  cubic  is  automatically  satisfied,  guaranteed  by  the  trigonometric  identity 
Eq.  A-53.  Therefore, 


39  =  cos  1  ^  4"  where  k  =  0,  ±1,  (A-60) 

and  from  Eqs.  A-54  and  A-57,  the  solutions  for  x  are 


Xk 


0,±1. 


(A-61) 


It  is  easy  enough  to  check  that  the  3  solutions  do  indeed  sum  to  zero  as  promised, 
since 

cos  9  +  cos 

for  arbitrary  9.  Two  of  the  solutions  will  be  positive,  corresponding  to  a  disk-like  and 
a  rod-like  cylinder,  and  the  third  solution  will  be  negative,  which  is  of  no  physical 
interest. 


„  2vr 

9  +  —  )  +  cos 
o 


(9-^'l=0 
3 


(A-62) 
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The  solutions  for  L/D  are 


/ !  o  / !  q 

—  =  .Xq  for  disk-like  cylinders  and  —  =  a:-,  1  for  rod-like  cylinders  in  Eq.  A-45; 


—  =  x0  3/2  for  disk-like  cylinders  and  —  =  aq  0/^  for  rod-like  cylinders  in  Eq.  A-46; 

I .  o  T .  o 

—  =  .x(y  for  disk- like  cylinders  and  —  =  for  rod-like  cylinders  in  Eq.  A-47. 


,-3/2  . 


A-3.2.  Cylinder  Probability  Distributions 


There  are  2  cases  to  consider:  rod- like  with  y  >  j  and  disk-like  with  y  <  y 


For  the  following  formulas,  define  a  =  (  — — 


-2/3 


7 T  L\  2//'?  71 


Dandb=UD)  4 


L/D  >  7t/4: 

The  PDF  is  given  by 


fii)  =  { 


q7  -  -  y2 

7max\/7max  7 
2  ay 


I  y2  /T2  —  0/2 

V  /max  V  /max  / 

and  the  CDF  is  given  by 


if  b  <  7  <  a 
if  a  <  7  <  7n 


f\  &7  +  «V7max  -  72  7  /  / 

1  1 - A; -  if  b  <  7  <  a 


^(7)  =  < 


72 

/max 


2a\/Tmax-72 


ry2 

I  max 


if  a  <  7  <  7„ 


(A-63) 


(A-64) 


Z/Z)  <  7r/4: 

The  PDF  is  given  by 

f  +  V 7max-72 


/  (7)  =  S 


7max\/7max-72 
2  ay 

I  y2  a/t2  —  y2 
V  /max  V  /max  / 
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if  a  <7  <6 
if  6  <  7  <  7n 


(A-65) 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 


and  the  CDF  is  given  by 


67  -  ay/ 7^ax  -  72 


^(7)  = 


72 

/  max 


1  - 


2a\/7 


2  —  /y2 

max  / 


These  functions  are  plotted  in  Fig.  A-4  for  the 


if  a  <  7  <  6 


if  b  <  7  <  7max 

case  when  L/D  =  1. 


(A-66) 


Fig.  A-4.  Plot  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  cylinder  with  L/D  —  1. 
Notice  the  discontinuity  at  7  =  (7t/4)-2/3  where  the  shape  factor  changes  from  being  single¬ 
valued  to  being  multivalued. 


Equations  A-64  and  A-66  can  be  solved  for  7  for  a  given  value  of  F — and  in 
this  way  turn  this  into  a  direct  and  fast  method  for  computing  the  shape  factor 
probability  distribution  of  a  random  tumbling  cylinder.  The  resulting  algorithm  has 
been  implemented  into  the  C++  code  in  Listing  A-3. 


Listing  A-3.  algo.cpp 

//  algo.cpp:  fast  algorithm  for  generating  the  shape  factor  of  an  RCC  with  a  uniform  random  orientation  over  the  unit 
sphere 

//  R.  Saucier,  Feb  2016  (see  p.  12  of  ARL-TR-2269  for  derivation  of  this  algorithm) 

#include  <iostream> 

#include  <cmath> 

#include  <chrono> 

#include  <random> 
using  namespace  std; 

int  main(  int  argc,  char*  argv[]  )  { 

double  l_d  =1.; 

if  (  argc  ==  2  )  l_d  =  atof(  argv[l]  );  //  L/D  =  1  is  default  or  override  with  1st  arg 

int  N  =  1000; 

if  (  argc  ==  3  )  N  =  atoi(  argv[2]  );  //  1000  samples  is  default  or  override  with  2nd  arg 

const  double  C  =  pow(  M_PI_4  *  l_d,  -2.  /  3.  ); 

const  double  A  =  C  *  l_d; 

const  double  B  =  C  *  M_PI_4; 

const  double  G_MAX2  =A*A+B*B; 

const  double  G_MAX  =  sqrt(  G_MAX2  ); 

const  double  K  =  G_MAX  /  (  2.  *  A  ); 

const  double  PI  =1.  -2.  *A*B/  G_MAX2; 
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24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 


const  double  P2 


=(B-A)*(B+A)/  G_MAX2 ; 


unsigned  seed  =  std: :chrono: : high_resolution_clock: :now() .time_since_epoch() .count() ; 
std : :mtl9937  rng(  seed  );  //  Mersenne  Twister  engine 

std : : uniform_real_distribution<double>  u(  0.,  1.  );  //  bind  uniform  distribution 

double  p,  q,  r,  g; 

for  (  int  i  =  0;  i  <  N;  i++  )  { 

p  =  u(  rng  ); 

q  =  1.  -  p; 

r  =  q  =4=  K; 

if  (  A  >=  B  )  {  //  same  as  L/D  >=  PI/4 

if  (  p<=Pl  )  g=B*q+A*  sqrt(  p  *  (  1.  +  q  )  ); 
else  g  =  G_MAX  *sqrt((l-r)*(l.  +  r)); 

> 

else  {  //  same  as  L/D  <  PI/4 

if  (  p<=P2  )  g=B*p+A*  sqrt(  q  *  (  1.  +  p  )  ); 
else  g  =  G_MAX  *  sqrt(  (  1.  -  r  )  *  (  1.  +r)); 

} 

std::cout  «  g  «  std::endl; 

} 

return  EXIT_SUCCESS; 


This  code  generates  over  28  million  shape  factors  per  second  on  a  Mac  with  a 
2.4-GHz  Intel  Xeon  processor.  Running  the  code  for  an  L/D  =  1  cylinder  gives  the 
results  shown  in  Fig.  A-5. 


Fig.  A-5.  Histograms  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  L/D  =  1  cylin¬ 
der  compared  to  plots  from  analytical  formulas,  Eqs.  A-65  and  A-66.  Notice  the  jump  in  the 
PDF  at  7  =  (7r/4)-2/3  sa  1.175,  as  predicted  (compare  to  Fig.  A-4). 
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A-4.  Tetrahedron 


Consider  a  regular  tetrahedron  of  unit  side  length,  which  has  a  volume  of  |  (^)3.  Its 
presented  area  ranges  from  to  The  PDF  as  a  function  of  area  is  given  by1 

12  .  ,/8T2-l\  8\/3  ,  /  3  —  20A2 \  8  1  „  1 

—  sin-1  - - —  +  -A—  cos-1  - —— - 7=  if  <A<  ^ 

7T  \l-4  A2  J  7 T  \3  —  16A2  J  Vs  Vs  V& 


f(A)=\6+™ 

Vs 


if  4=  <A<^ 
\/6  “  "4 

if^<T<i 

4  ~  ~  2 


(A-67) 


Since  presented  area,  A,  and  dimensionless  shape  factor,  7,  are  related  by 


A  =  7C2/3, 


(A-68) 


where  V  is  the  volume,  we  can  get  the  PDF  as  a  function  of  7  by  simply  scaling  the 
area  by  V~2/3  =  2  ■  32/3. 

The  CDF  as  a  function  of  area  is  given  by 

12  „  .  ,/8T2-l\  8^3,  ,  / 3  —  20T2 \  8  „ 

— A  sm_1  - —  H - Tcos_i  - — - =A+ 

7 r  \1-4A2)  7 r  \3  -  16A2  J  VS 


—  tan 
7 r 


—  tan 
7 r 


V2(l  -  3 A) 
Vl  -  6A2 


V~2.(2  -  3 VSA) 
Vl  -  6 A2 


_!  f  v/2(l  +  3T) 
\  Cl  -  6A2 


C2(2  +  3C3A) 
Cl  -  6A2 


J6  r 

—  tan_1(3  —  2C2)  —  tan^1(3  +  2C2)  — 

7T  L 

—  tan_1(4C2  —  3C3)  +  tan_1(4C2  +  3C3) 

7T  L 

^6+-^0A+  —  tan_1(3  —  2C2)  —  tan_1(3  +  2C2) 
—  tan_1(4C2  —  3C3)  +  tan_1(4C2  +  3C3) 

7T  L 


6T-2 


if  ~^<A<^= 
C8  C6 


VB  -  4 

4  ~  ~  2 

(A-69) 


The  C++  code  in  Listing  A-4  implements  the  PDF  and  CDF,  Eqs.  A-67  and  A-69. 
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Listing  A-4.  tetrahedron.cpp 


1  //  tetrahedron.cpp:  generates  plotting  points  for  tetrahedron  pdf  and  cdf 

2  //  Ref:  Vickers,  G.  T.  and  Brown,  D.  J., 

3  //  "The  distribution  of  projected  area  and  perimeter  of  convex,  solid  particles," 

4  //  Proc.  R.  Soc.  Lond.  A  (2001),  Vol.  457,  pp.  283-306. 

5 

6  #include  <iostream> 

7  #include  <cmath> 

8  #include  <cstdlib> 

9  #include  <cassert> 

10 

11  const  double  CUBE- VOLUME  =  1.  /  (  M-SQRT2  *  M-SQRT2  *  M-SQRT2  ); 

with  unit  side  length 

12  const  double  TETRAHEDRON- VOLUME  =  CUBE- VOLUME  /  3.; 

13  const  double  FACTOR  =  pow(  TETRAHEDRON. VOLUME ,  -2.  /  3.  ); 

dimensionless  shape  factor 

14  const  double  SQRT3  =  sqrt(  3.  ); 

15  const  double  SQRT6  =  M.SQRT2  *  SQRT3 ; 

16  const  double  SQRT8  =  sqrt(  8.  ); 

17  const  double  AMIN  =  1.  /  SQRT8 ; 

18  const  double  AMAX  =  1.  /  2.; 

19 

20  double  pdf(  double  x  )  { 

21 

22  assert (  AMIN  <=  x  &&  x  <=  AMAX  ); 

23 

24  if  (  x  <  1.  /  SQRT6  )  { 

25  double  x2  =  x  *  x; 

26  return  (  12.  /  M_PI  )  *  asin(  (  8.  *  x2 

27  (  8.  *  SQRT3  /  M_PI  )  *  acos(  (  3 

28  > 

29  else  if  (  x  <  SQRT3  /  4.  ) 

30  return  6.  +  16.  /  SQRT3 ; 

31  else 

32  return  6. ; 

33  > 

34 

35  double  cdf(  double  x  )  { 

36 

37  assert (  AMIN  <=  x  &&  x  <=  AMAX  ); 

38 

39  if  (  x  <  1.  /  SQRT6  )  { 

40 

41  double  x2  =  x  *  x; 

42  double  x4  =  x2  *  x2; 

43  return 


44 

(  -2. 

*  (  4. 

*  SQRT3  *  M_PI 

*  x  -  12.  *  SQRT3  *  x  *  acos(  (3.  -  20.  *  x2  )  /  (  3.  -  16.  *  x2  )  )+ 

45 

18. 

*  x  * 

asin(  (1.  -  8. 

*  x2  )  /  (  1.  -  4.  *  x2  )  )  - 

46 

9. 

* 

atan( 

3.  -  2.  *  M-SQRT2  )  + 

47 

9. 

* 

atan( 

3.  +2.  *  M.SQRT2  )  + 

48 

9. 

* 

atan( 

4.  *  M.SQRT2  - 

3.  *  SQRT3  )  + 

49 

9. 

* 

atan( 

4.  *  M.SQRT2  + 

3.  *  SQRT3  )  - 

50 

9. 

* 

atan( 

(  M_SQRT2  *  ( 

1.  -  3.  *  x  )  )  /  sqrt(  1.  -  6.  *  x2  )  )  - 
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9. 

* 

atan( 

(  M_SQRT2  *  ( 

1.  +  3.  *  x  )  )  /  sqrt(  1.  -  6.  *  x2  )  )  - 
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9. 

* 

atan( 

(  M-SQRT2  *  ( 

-2.  +  3.  *  SQRT3  *  x  )  *  sqrt(  x2  -  6.  *  x4  )  )  /  (  x  *  (  -1.  +  6.  *  x2  )  ) 

)  + 
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9. 

* 

atan( 

(  M-SQRT2  *  ( 

2.  +3.  *  SQRT3  *  x  )  *  sqrt(  x2  -  6.  *  x4  )  )/(  x  *  (  -1.  +  6.  *  x2  )  )  ) 

)  )  /  (  3. 

M_PI 

); 
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> 

55 

else  if 

( 

x  <  SQRT3  /  4.  )  { 

56 

return 

(  6. 

+  16.  /  SQRT3  ) 

*  x  + 

57  (  6.  *  (  atan(  3.  -  2.  *  M.SQRT2  )  - 

58  atan(  3.  +2.  *  M.SQRT2  )  - 

59  atan(  4.  *  M.SQRT2  -  3.  *  SQRT3  )  - 

60  atan(  4.  *  M.SQRT2  +  3.  *  SQRT3  )  )  )  /  M_PI; 

61  } 

62  else 

63  return  6.  *  x  -  2. ; 

64  > 

65 

66  int  main(  int  argc,  char**  argv  )  {  //  specify  number  of  points  on  commandline  or  use  1000 

67 

68  int  N  =  1000; 

69  if  (  argc  ==  2  )  N  =  atoi(  argv[l]  ); 

70 

71  for  (  double  a  =  AMIN;  a  <=  AMAX;  a  +=  (  AMAX  -  AMIN  )  /  double (  N  )  )  { 

72 

73  double  g  =  FACTOR  *  a; 

74  std::cout  «  g  «  "\t" 

75  «  pdf (  a  )  /  FACTOR  «  "\t" 

76  «  cdf(  a  )  «  std::endl; 

77  > 

78  return  EXIT-SUCCESS; 

79  > 


-  1.  )  /  (  1.  -  4.  *  x2  )  )  + 

.  -  20.  *  x2  )  /  (  3.  -  16.  *  x2  )  )  -  8.  /  SQRT3 ; 


//  volume  of  cube  that  encloses  the  tetrahedron 

//  volume  of  tetrahedron  with  unit  side  length 
//  factor  that  converts  presented  area  to 
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This  is  plotted  in  Fig.  A-6. 


Fig.  A-6.  Plot  of  shape  factor  PDF  and  CDF  for  a  randomly  oriented  regular  tetrahedron 


A-5.  Ellipsoid 

The  equation  of  an  ellipsoid,  in  which  its  principal  axes  are  aligned  with  the  coordi¬ 
nate  axes,  is  given  by 


x 2  y 2  z 2 

- h  —  H - =  1 

a 2  b2  c2 


(A-70) 


where  a,  b,  and  c  are  the  3  semi -principal  axes,  and  its  volume  is 


V  =  -nabc. 
3 


(A-71) 


Let  the  viewing  angle  be  specified  by  the  polar  angle  9  and  the  azimuthal  angle  0. 
The  presented  area  of  the  ellipsoid  will  be  in  the  shape  of  an  ellipse  of  area5 


A  =  n 


sin2  9  cos2  0  +  a2c2  sin2  6  sin2  0  +  a2b 2  cos2  9. 


(A-72) 


Let  the  dimensions  be  ordered  so  that  a  <  b  <  c.  The  CDF  as  a  function  of  area  is 
given  by1 


F(x)  =  {  l 

1  7T 


i  -  -Vin 


Rf{ 0,£  -  !,£-»?)-  g-Rj(0,£  -  1,£  -  r/,0 


Sin  1  \/(A2m  ~  ^?nin)/(Aiax  -  ^7) 


i  -  -VCn 


Rf{ 0, 1  -  C 1  -  v)  -  3^j(0, 1  -  C 1  -  y,  i) 


if  dmiii  7  ^  A  I, 

if  x  =  Am 


if  Am  <  x  <  An 


(A-73) 


5Vickers  GT.  The  projected  areas  of  ellipsoids  and  cylinders.  Powder  Technology.  1996;  86: 
195-200. 
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where 


£  -  ;rX  A2  ’  ??=  ^"aX  A2  .  4nin  =  7T0&,  Am  =  7TOC,  Amax  =  nbc,  (A-74) 

^max  ^min 


max  m 


and  the  functions  7iV  and  f?j  are  the  Carlson  symmetrized  form  of  the  classic  elliptic 
integrals: 

i  r°°  dt 

(A-75) 


„ ,  x  i  r  dt 

nF(x,y,z)  =  -  /  — .  = 

2  Jo  y/(t  +  x)(t  +  y)(t  +  z) 


and 


3  r°° 

Rj(x,y,z,p)  =  - 


dt 


(A-76) 


2  7o  (7  +  p)  a/(7  +  x)(7  +  p)(7  +  z) 

It  is  shown  in  Numerical  Recipes 6  that  the  PDF  as  a  function  of  area  is  given  by 


/  0)  = 


2x 

Rf(0,  (Am  —  x  )(Amax  —  Amin),  (Amax  —  x  )(Am  —  Amin))  if  Amjn  <  x  <  Am 
7 r 

2a; 

—Rf( 0,  (a:  —  Aia)(Amax  —  Amin),  (a;  —  Amin)(Amax  —  Am))  if  Am  <  x  <  Amax 

(A-77) 

There  is  a  logarithmic  singularity  in  /  at  x  =  Am  when  a,  6,  and  c  are  all  different. 
The  total  surface  area  of  the  ellipsoid  is  given  by 

„  „  2  2t rbc  r,  (  1  1  1  \  2i:abc  (  1  1  \  (  1  1  \  „  /  1  1  1  1 


~2  .2  )  l  „2  b2  )  RJ  [  c2  »  b2  >  >  fl2  )  ’ 


(A-78) 

and  from  this  we  can  get  the  mean  presented  area  as  5/4  from  Cauchy’s  theorem. 

The  presented  areas  in  these  formulas  are  easily  converted  to  dimensionless  shape 
factors,  7,  via  the  relationship 

A  =  jV2/3,  (A-79) 

where  V  is  the  ellipsoid  volume.  To  convert  the  PDF  from  a  function  of  area  to  a 
function  of  shape  factor,  use 


r  /  \  dF  dF  d A  -tr2/3r/  \ 

/h)  =  77  =  iI77  =  l/  /W' 


(A-80) 


The  formulas  in  Eqs.  A-77  and  A-74  are  plotted  in  Fig.  A-7. 


6Press  WH,  Teukolsky  SA,  Vetterling  WT,  Flannery  BP.  Numerical  recipes  in  C:  the  art  of 
scientific  computing.  New  York  (NY):  Cambridge  University  Press;  1995. 
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Fig.  A-7.  Plot  of  shape  factor  PDF  and  CDF  for  an  ellipsoid  with  a  =  1,  b  =  2,  c  =  3.  There 
is  a  logarithmic  singularity  in  the  PDF  when  a,  h,  and  c  are  all  different  and  is  located  at 

7  =  (4/3) _2//3(7rac/h2)1/3,  which  in  this  case  is  at  37r1//3/4  «  1.098. 

Listings  A-5  through  A-8  can  be  used  for  printing  out  the  solid  curves  in  Fig.  A-7. 


Listing  A-5.  ellipsoid.cpp 

//  ellipsoid.cpp:  Prints  out  the  PDF  and  CDF  as  a  function  of  shape  factor  for  an  ellipsoid 
//  Specify  the  dimension  of  the  ellipsoid  on  the  commandline  (or  use  default  3/2/1) 

//  R.  Saucier,  Feb  2016 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

double  rf(  double  x,  double  y,  double  z  ); 

double  r j (  double  x,  double  y,  double  z,  double  p  ); 

inline  double  min(  double  a,  double  b,  double  c  )  {  return  std::min(  std::min(  a,  b  ),  c  );  } 

inline  double  max(  double  a,  double  b,  double  c  )  {  return  std::max(  std::max(  a,  b  ),  c  );  } 

inline  double  mid(  double  a,  double  b,  double  c  )  {  return  std::max(  std::min(  a,  b  ),  std::min(  std::max(  a,  b  ),  c  )  ); 

} 

class  PDF  {  //  functor  for  probability  density  function 

public: 

PDF(  double  amin,  double  am,  double  amax  )  :  _amin(  amin  ),  _am(  am  ),  _amax(  amax  )  {> 
double  operator()(  double  x  )  { 
double  y,  z; 

if  (  _amin  <=  x  &&  x  <=  _am  )  { 

y  =  (  _am  -  x  )  *  (  _am  +  x  )  *  (  _amax  -  _amin  )  *  (  _amax  +  _amin  ) ; 

z  =  (  _amax  -  x  )  *  (  _amax  +  x  )  *  (  _am  -  _amin  )  *  (  _am  +  _amin  ) ; 

} 

else  { 

y  =  (  x  -  _am  )  *  (  x  +  _am  )  *  (  _amax  -  _amin  )  *  (  _amax  +  _amin  ) ; 

z  =  (  x  -  _amin  )  *  (  x  +  _amin  )  *  (  _amax  -  _am  )  *  (  _amax  +  _am  ) ; 

} 

return  (  x  /  M_PI_2  )  *  rf(  0,  y,  z  ); 

> 

private: 

double  _amin,  _am,  _amax; 

}; 


class  CDF  {  //  functor  for 

cumulative  distribution  function 

public: 

CDF(  double  amin,  double 

am,  double  amax  )  :  _amin(  amin  ) 

,  _am(  am  ),  _amax(  amax  ) 

{> 

double  operator()(  double 

x  )  { 

if  (  x  ==  _am  ) 

return  asin(  sqrt( 

(  _am  -  _amin  )  *  (  _am  +  _amin  ) 

/  (  (  _amax  -  _amin  )  *  ( 

'  _amax  +  _amin  )  ) 

)  )  /  M_PI_2 ; 

double 

y,  z,  p,  xi,  eta; 

if  (  _i 

amin  <=  x  &&  x  <=  _am  )  { 

xi 

=  (  _amax  -  x  )  *  (  _amax 

+ 

x  : 

i  /  ( 

(  _amax  - 

_am  ; 

1  *  ( 

_amax 

+  _am  )  ) ; 

eta 

=  (  _amax  -  x  )  *  (  _amax 

+ 

x  : 

i  /  ( 

(  _amax  - 

_amin  ; 

1  *  ( 

_amax 

+  _amin  )  ) ; 

y 

=  xi  -  1.; 

z 

=  xi  -  eta; 

P 

=  xi; 

} 

else  { 

xi 

=  (  _amax  -  x  )  *  (  _amax 

+ 

x  : 

i  /  ( 

(  _amax  - 

_am  ; 

1  *  ( 

_amax 

+  _am  )  ) ; 

eta 

=  (  _amax  -  x  )  *  (  _amax 

+ 

x  ] 

i  /  ( 

(  _amax  - 

_amin  ; 

1  *  ( 

_amax 

+  _amin  )  ) ; 

y 

=  1.  -  xi; 

z  =  1.  -  eta; 
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81 
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84 
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86 
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P  =1; 

} 

return  1.  -  (  sqrt(  xi  *  eta  )  /  M_PI_2  )  *  (  rf(  0,  y,  z  )  -  rj(  0,  y,  z,  p  )  /  3.  ); 

> 

private: 

double  _amin,  _am,  _amax; 

}; 


int  main(  int  argc,  char*  argv[]  )  { 


double  a  =  1 

,  b  =  2.  ,  c  =  3. ; 

if  (  argc  == 

4  )  {  //or  specify  the  3 

dimensions  (in  any  order)  on  the  command  line 

a  =  atof( 

argv[l]  ); 

b  =  atof( 

argv[2]  ); 

c  =  atof( 

argv[3]  ); 

} 

const  double 

A  =  min(  a,  b,  c  ) ; 

//  minimum  value 

const  double 

B  =  mid(  a,  b,  c  ) ; 

//  intermediate  value 

const  double 

C  =  max(  a,  b,  c  ) ; 

//  maximum  value 

const  double 

V  =  (  4.  /  3.  )  *  M_PI  *  A  * 

B  *  C;  //  ellipsoid  volume 

const  double 

SI  =  pow (  V,  +2.  /  3.  ); 

//  factor  to  convert  PDF  from  area  to  shape  factor 

const  double 

S2  =  pow (  V,  -2.  /  3.  ); 

//  factor  to  convert  area  to  shape  factor 

const 

double 

AMIN  =  M_PI 

* 

A 

* 

B; 

const 

double 

AM  =  M_PI 

* 

A 

* 

C; 

const 

double 

AMAX  =  M_PI 

* 

B 

* 

C; 

PDF  pdf (  AMIN,  AM,  AMAX  ); 

CDF  cdf (  AMIN,  AM,  AMAX  ); 
const  int  N  =  1000; 

for  (  double  x  =  AMIN;  x  <=  AMAX;  x  +=  (  AMAX  -  AMIN  )  /  double!  N  )  ) 

std::cout  «  S2  *  x  «  "\t"  «  SI  *  pdf (  x  )  «  "\t"  «  cdf(  x  )  «  std::endl; 

return  EXIT.SUCCESS; 


Listing  A-6.  rf.cpp 


//  rf.cpp: 

Computes 

Carlson's  elliptic  integral  of  the  first 

kind,  Rf(x,y,z), 

// 

where  x, 

y,  and  z  must  be  nonnegative  and  at  most 

one  can  be  zero.  TINY 

// 

must  be 

at  least  5  times  the  machine  underflow  limit  and  BIG  at  most 

// 

one  fifth  the  machine  overflow  limit. 

//  Ref:  Press,  W.H.,  Teukolsky,  S.A.,  Vetterling,  W,  T.,  Flannery,  B.P., 

//  Numerical  Recipes  in  C,  Cambridge  University  Press,  1992. 

#include  <cmath> 

#include  <cstdlib> 

#include  <iostream> 

inline  double  FMIN3(  double  a,  double  b,  double  c  )  {  return  std::min(  std::min(  a,  b  ),  c  );  } 

inline  double  FMAX3(  double  a,  double  b,  double  c  )  {  return  std::max(  std::max(  a,  b  ),  c  );  } 

double  rf(  double  x,  double  y,  double  z  )  { 

const  double  ERRTOL  =  0.08; 

const  double  TINY  =  5.7e-103;  //  at  least  2(DBL_MIN)/'(l/3) ,  where  DBL-MIN  =  2.22507e-308 

const  double  BIG  =  l.le+102;  //  at  most  (1/5) (DBL.MAXJ'Ml/S) ,  where  DBL_MAX  =  1.79769e+308 

const  double  THIRD  =  1.  /  3.; 

const  double  Cl  =  1.  /  24.; 

const  double  C2  =0.1; 

const  double  C3  =  3.  /  44.; 

const  double  C4  =  1.  /  14.; 

double  alamb,  ave,  delx,  dely,  delz,  e2,  e3,  sqrtx,  sqrty,  sqrtz,  xt,  yt,  zt; 


if  (  FMIN3 (  x,  y,  z  )  <  0.0  |  | 

FMIN3 (  x+y,  x+z,  y+z)<  TINY  | | 
FMAX3 (  x,  y,  z  )  >  BIG  )  { 


std : 

cerr  « 

"invalid 

arguments  in 

rf :  "  «  std : : endl ; 

std : 

cerr  « 

"  x  = 

«  x  «  std 

:  endl 

« 

"  y  = 

«  y  «  std 

:  endl 

« 

"  z  = 

«  z  «  std 

: endl ; 

exit(  EXIT. FAILURE  ); 

} 


xt  =  x; 
yt  =  y; 
zt  =  z; 
do  { 


sqrtx  =  sqrt(  xt  ) 

sqrty  =  sqrt(  yt  ) 

sqrtz  =  sqrt(  zt  ) 

alamb  =  sqrtx  *  (  sqrty  +  sqrtz  )  +  sqrty  *  sqrtz; 
xt  =  0.25  *  (  xt  +  alamb  ); 
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47 

yt  =  0.25  *  (  yt  +  alamb  ); 

48 

zt  =  0.25  *  (  zt  +  alamb  ); 

49 

ave  =  THIRD  *  (  xt  +  yt  +  zt  ) ; 

50 

delx  =  (  ave  -  xt  )  /  ave; 

51 

dely  =  (  ave  -  yt  )  /  ave; 

52 

delz  =  (  ave  -  zt  )  /  ave; 

53 

}  while (  FMAX3 (  fabs(  delx  ),  fabs( 

dely  ),  fabs(  delz  ) 

)  >  ERRTOL  ); 

54 

e2  =  delx  *  dely  -  delz  *  delz; 

55 

e3  =  delx  *  dely  *  delz; 

56 

return  (  1.  +  (  Cl  *  e2  -  C2  -  C3  * 

e3  )  *  e2  +  C4  *  e3  ] 

\  /  sqrt(  ave  ); 

57 

> 

Listing  A-7.  rj.cpp 


1  //  rj.cpp:  Computes  Carlson's  elliptic  integral  of  the  third  kind,  Rj(x,y,z,p), 

2  //  where  x,  y,  and  z  must  be  nonnegative  and  at  most  one  can  be  zero,  p 

3  //  must  be  nonzero.  If  p  <  0,  the  Cauchy  principal  value  is  returned.  TINY 

4  //  must  be  at  least  twice  the  cube  root  of  the  machine  underflow  limit  and 

5  //  BIG  at  most  one  fifth  the  cube  root  of  the  machine  overflow  limit. 

6  //  Ref:  Press,  W.H.,  Teukolsky,  S.A.,  Vetterling,  W,  T.,  Flannery,  B.P., 

7  //  Numerical  Recipes  in  C,  Cambridge  University  Press,  1992. 

8 

9  #include  <cmath> 

10  #include  <cstdlib> 

11  #include  <iostream> 

12 

13  double  rc(  double  x,  double  y  ); 

14  double  rf(  double  x,  double  y,  double  z  ); 

15 

16  inline  double  FMIN3(  double  a,  double  b,  double  c  )  {  return  std::min(  std::min(  a,  b  ),  c  );  } 

17  inline  double  FMAX3(  double  a,  double  b,  double  c  )  {  return  std::max(  std::max(  a,  b  ),  c  );  } 

18  inline  double  FMIN4(  double  a,  double  b,  double  c,  double  d  )  {  return  std::min(  FMIN3(  a,  b,  c  ),  d  );  } 

19  inline  double  FMAX4(  double  a,  double  b,  double  c,  double  d  )  {  return  std::max(  FMAX3(  a,  b,  c  ),  d  );  } 

20  inline  double  SQR(  double  a  )  {  return  a  *  a;  } 

21 

22  double  r j (  double  x,  double  y,  double  z,  double  p  )  { 

23 

24  const  double  ERRTOL  =  0.05; 


25 

const 

double 

TINY 

= 

5.7e-103; 

// 

at 

least 

2(DBL_MINr(l/3) ,  where  DBL.MIN  =  2.22507e-308 

26 

const 

double 

BIG 

= 

l.le+102; 

// 

at 

most 

(1/5)  (DBL_MAX)yv(l/3) ,  where  DBL_MAX  =  1.79769e+308 

27 

const 

double 

Cl 

= 

3.  /  14.; 

28 

const 

double 

C2 

= 

1.  /  3.; 

29 

const 

double 

C3 

= 

3.  /  22.; 

30 

const 

double 

C4 

= 

3.  /  26.; 

31 

const 

double 

C5 

= 

0.75  *  C3 ; 

32 

const 

double 

C6 

= 

1.5  *  C4; 

33 

const 

double 

C7 

= 

0.5  *  C2; 

34 

const 

double 

C8 

= 

C3  +  C3; 

35 

36  double  a  =  0.,  alamb,  alpha,  ans,  ave,  b  =  0.,  beta,  delp,  delx,  dely,  delz,  ea,  eb,  ec, 

37  ed,  ee,  fac,  pt,  rex  =0.,  rho,  sqrtx,  sqrty,  sqrtz,  sum,  tau,  xt,  yt,  zt; 

38 

39  if  (  FMIN3 (  x,  y,  z  )  <  0.0  |  | 

40  FMIN4(  x+y,  x+z,  y+z,  fabs(  p  )  )  <  TINY  | | 

41  FMAX4(  x,  y,  z,  fabs(  p  )  )  >  BIG  )  { 


42 

std : : cerr  «  " 

invalid 

arguments  in  rj :  "  «  std::endl; 

43 

std : : cerr  «  " 

X 

=  "  «  x  «  std : : endl 

44 

«  11 

y 

=  11  «  y  «  std : :  endl 

45 

«  " 

z 

=  "  «  z  «  std : : endl 

46 

«  " 

p 

=  "  «  p  «  std : : endl 

47 

«  " 

TINY 

=  "  «  TINY  «  std : :  endl 

48 

«  " 

BIG 

=  11  «  BIG  «  std:: endl; 

49 

exit(  EXIT. FAILURE  ); 

50 

} 

51 

52 

sum  =  0. ; 

53 

fac  =  1. ; 

54 

if  (  P  >  0.  )  { 

55 

xt  =  x; 

56 

yt  =  y; 

57 

zt  =  z; 

58 

pt  =  p; 

59 

} 

60 

else  { 

61 

xt  =  FMIN3 (  x, 

y,  z  ) 

62 

zt  =  FMAX3 (  x, 

y,  z  ) 

63 

yt  =  x  +  y  +  z 

-  xt  - 

zt; 

64 

a  =  1.  /  (  yt 

-  p  ); 

65 

b  =  a  *  (  zt  - 

yt  )  * 

(  yt  -  xt  ); 

66 

pt  =  yt  +  b; 

67 

rho  =  xt  *  zt 

/  yt; 

68 

tau  =  p  *  pt  / 

yt; 

69 

rex  =  rc(  rho, 

tau  ) ; 
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27 

28 
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} 


do  { 

sqrtx  =  sqrt(  xt  ) 

sqrty  =  sqrt(  yt  ) 

sqrtz  =  sqrt(  zt  ) 

alamb  =  sqrtx  *  (  sqrty  +  sqrtz  )  +  sqrty  *  sqrtz; 

alpha  =  SQR(  pt  *  (  sqrtx  +  sqrty  +  sqrtz  )  +  sqrtx  *  sqrty  *  sqrtz  ); 
beta  =  pt  *  SQR(  pt  +  alamb  ); 
sum  +=  fac  *  rc(  alpha,  beta  ); 
fac  *=  0.25; 


xt  =  0.25  *  ( 

xt  + 

alamb  ) ; 

yt  =  0.25  *  ( 

yt  + 

alamb  ) ; 

zt  =  0.25  *  ( 

zt  + 

alamb  ) ; 

pt  =  0.25  *  ( 

pt  + 

alamb  ) ; 

ave  =  0.2  *  ( 

xt  + 

yt  +  zt  +  pt  +  pt  ) ; 

delx  =  (  ave 

xt 

/  ave; 

dely  =  (  ave 

yt 

/  ave; 

delz  =  (  ave 

zt 

/  ave; 

delp  =  (  ave 

Pt 

/  ave; 

}  while (  FMAX4(  fabs( 

delx  ),  fabs(  dely  ),  fabs(  delz  ),  fabs(  delp  ) 

)  >  ERRTOL  ); 

ea  =  delx  *  (  dely  + 

delz  )  +  dely  *  delz; 

eb  =  delx  *  dely 

*  delz; 

ec  =  delp  *  delp; 

ed  =  ea  -  3.  *  ec; 

ee  =  eb  +  2.  *  delp  * 

(  ea  -  ec  ) ; 

ans  =  3.  *  sum  + 

fac  *  (  1. 

+  ed 

*  (  -Cl  +  C5  *  ed  -  C6  *  ee  )  + 

eb 

*  (  C7  +  delp  *  (  -C8  +  delp  *  C4  )  )  + 

delp  *  ea  *  (  C2  -  delp  *  C3  )  - 

C2 

*  delp  *  ec  )  /  (  ave  *  sqrt(  ave  )  ); 

if  (  p  <=  0.  )  ans  = 

*  (  b  *  ans  +  3.  *  (  rex  -  rf(  xt,  yt,  zt  )  )  ) 

return  ans; 

} 


Listing  A-8.  rc.cpp 

//  rc.cpp:  Computes  Carlson's  degenerate  elliptic  integral,  Rc(x,y),  where  x  must 
//  be  nonnegative  and  y  must  be  nonzero.  If  y  <  0,  the  Cauchy  principal 

//  value  is  returned.  TINY  must  be  at  least  5  times  the  machine  underflow 

//  limit  and  BIG  must  be  at  most  one  fifth  the  machine  overflow  limit. 

//  Ref:  Press,  W.H.,  Teukolsky,  S.A.,  Vetterling,  W,  T.,  Flannery,  B.P., 

//  Numerical  Recipes  in  C,  Cambridge  University  Press,  1992. 


#include  <cmath> 
#include  <cstdlib> 
#include  <iostream> 


double  rc(  double  x,  double  y  )  { 
const  double  ERRTOL  =  0.04; 

const  double  TINY  =  5.7e-103;  //  at  least  2(DBL_MIN)/'(l/3) ,  where  DBL_MIN  =  2.22507e-308 

const  double  BIG  =  l.le+102;  //  at  most  (1/5) (DBL_MAX)/'(l/3) ,  where  DBL_MAX  =  1.79769e+308 

const  double  SQRTNY  =  sqrt(  TINY  ); 
const  double  TNBG  =  TINY  *  BIG; 

=  2.236  /  SQRTNY; 

=  TNBG  *  TNBG  /  25.; 

=  1.  /  3.; 

=  0.3; 

=  1.  /  7.; 

=  0.375; 

=  9.  /  22.; 


const  double  C0MP1 
const  double  C0MP2 
const  double  THIRD 
const  double  Cl 
const  double  C2 
const  double  C3 
const  double  C4 


double  alamb,  ave,  s,  w,  xt,  yt; 
if  (  x  <  0.0  ||  y  ==  0.0  || 

(  x  +  fabs (  y  )  )  <  TINY  | | 

(  x  +  fabs (  y  )  )  >  BIG  | I 
(  y  <  -C0MP1  &&x>0.0&&x<  C0MP2  )  )  { 
std::cerr  «  "invalid  arguments  in  rc:  "  «  std::endl 
«  "  x  =  "  «  x  «  std : :  endl 

«  "  y  =  "  «  y  «  std : : endl ; 

exit(  EXIT_ FAILURE  ); 

> 

if  (  y  >  0.  )  { 
xt  =  x; 
yt  =  y; 
w  =  1. ; 

} 

else  { 

xt  =  x  -  y; 
yt  =  -y; 

w  =  sqrt(  x  )  /  sqrt(  xt  ); 

} 

do  { 
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48 

alamb  =  2.  *  sqrt( 

xt  )  *  sqrt(  yt  )  +  yt; 

49 

xt  =  0.25  *  (  xt  + 

alamb  ) ; 

50 

yt  =  0.25  *  (  yt  + 

alamb  ) ; 

51 

ave  =  THIRD  *  (  xt 

+  yt  +  yt  ); 

52 

s  =  (  yt  -  ave  )  / 

ave; 

53 

}  while  (  fabs(  s  )  > 

ERRTOL  ); 

54 

return  w  *  (  1.  +  s  * 

s*(Cl+s*(C2+s*i 

(  C3  +  s  *  C4  )  ) 

)  )  /  sqrt (  ave  ); 

55 

} 
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Appendix  B.  Uniform  Sampling  over  the  Unit  Sphere 
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The  global  version  of  Archimedes’  theorem1  states  that  the  area  of  a  sphere  is  equal 
to  the  area  of  a  cylinder  circumscribed  about  the  sphere,  excluding  the  bases.  The  area 
of  a  unit  sphere  is  47t.  The  area  of  the  circumscribed  cylinder  is  the  circumference 
times  the  height:  27t  x  2  =  Ait.  The  local  version  of  the  theorem  states  further  that 
any  region  on  the  sphere  is  equal  to  the  axial  projection  on  the  cylinder.  This  is  a 
very  powerful  theorem  for  our  purposes  since  it  is  much  easier  to  define  a  sampling 
strategy  on  the  cylinder,  which  we  can  lay  out  flat  and  independently  sample  0  and 
z,  and  then  use  Archimedes’  theorem  to  map  onto  the  unit  sphere. 

Let  6  and  0  be  the  polar  and  azimuthal  angles,  respectively,  on  the  unit  sphere,  and  let 
0  and  z  be  coordinates  on  the  circumscribed  cylinder,  where  9  G  [0,  n],  0  G  [0,  27t], 
and  z  G  [—1,1].  Then  the  mapping  from  the  cylinder  to  the  sphere  [0,  27r]  x  NM] 
=>•  S2(9,  0)  is  simply 

9  =  cos-1 z  (B-l) 

while  the  0  value  remains  the  same.  Now  that  we  know  the  mapping  from  the 
cylinder  to  the  sphere,  we  focus  on  the  sampling  strategy  on  the  unwrapped  cylinder, 
the  aspect  ratio  of  which  is  depicted  in  Fig.  B-l. 

+1 - 

z 

-1 - 

0  0  27T 

Fig.  B-l.  Sampling  on  the  [0,  2-n]  x  [—1,1]  circumscribed  cylinder  allows  us  to  sample  both 
0  and  z  uniformly  and  independently  over  their  entire  range 


We  describe  4  sampling  strategies,  2  randomized  and  2  deterministic. 


'Shao  M,  Badler  N.  Spherical  sampling  by  Archimedes’  theorem.  Philadelphia  (PA):  University 
of  Pennsylvania:  1996;  Technical  Report  MS-CIS-96-02. 
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B-1.  Uniform  Random 


The  first  is  independent  uniform  random  sampling2  on  both  0  ~  [7(0,  27t)  and 
z  ~  U{— 1, 1).  The  C++  code  is  given  in  Listing  B-1. 

Listing  B-1.  uniform.cpp 


//  uniform.cpp:  generate  a  uniform  random  distribution  over  the  unit  sphere 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

#include  <iomanip> 

int  main(  int  argc,  char*  argv[]  )  { 

unsigned  int  N  =  1000; 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  ); 

std::cout  «  std : : setprecision (6)  «  std:: fixed; 

rng : : Random  rng; 

double  th,  ph,  x,  y,  z; 

for  (  unsigned  int  n  =  0;  n  <  N;  n++  )  { 

x  =  rng.uniform(  0.,  2.  *  M_PI  ); 
z  =  rng.uniform(  -1.,  1  ); 

//std::cout  «  x  «  "\t"  «  z  «  std::endl; 
ph  =  x; 

th  =  acos(  z  ) ; 

//std::cout  «  ph  «  "\t"  «  th  «  std::endl; 
x  =  sin(  th  )  *  cos(  ph  ) ; 
y  =  sin(  th  )  *  sin(  ph  ) ; 

std::cout  «  x  «  "\t"  «  y  «  "\t"  «  z  «  std::endl; 

> 

return  EXIT_SUCCESS; 


This  is  the  simplest  strategy,  and  it  has  the  advantage  that  we  do  not  need  to  know 
the  total  number  of  sample  points  beforehand;  we  can  simply  continue  until  we  meet 
some  convergence  criterion.  The  biggest  disadvantage  is  that  it  produces  a  pattern 
that  contains  clustering  of  points  and  is  not  very  uniform,  as  we  see  in  Fig.  B-2. 

B-2.  Stratified  Random 

We  can  improve  upon  the  clustering  problem  that  we  see  with  uniform  sampling  by 
using  stratified  random  sampling.  This  can  be  achieved  by  imposing  a  grid  on  the 
cylinder  and  drawing  a  random  sample  within  each  grid  cell.  The  C++  code  is  given 
in  Listing  B-2. 


Listing  B-2.  strat.cpp 


//  strat.cpp:  Stratified  uniform  spherical  sampling  based  upon  Archimedes'  Theorem; 

//  works  by  subdividing  the  [0,2  pi]  x  [-1,1]  cylinder  into  N~2  rectangles, 

//  randomly  selects  a  point  from  each,  and  then  maps  onto  the  unit  sphere. 

//  Reference:  Min-Zhi  Shao  and  Norman  Badler,  "Spherical  Sampling  by  Archimedes'  Theorem, 
//  http://repository.upenn.edu/cis_reports/184/,  25  June  2007. 

#include  "Random. h" 


2Saucier  R.  Computer  generation  of  statistical  distributions.  Aberdeen  Proving  Ground  (MD): 
Army  Research  Laboratory  (US);  2000  Mar.  Report  No.:  ARL-TR-2168. 
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#include  <iostream> 
#include  <cstdlit» 
#include  <cmath> 


int  main(  int  argc,  char*  argv[]  )  {  //  override  default  N  on  commandline 

int  N  =  32;  //  number  of  points  is  N/'2,  so  default  is  32^2  =  1024 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  ); 

const  double  DEL_X  =  2.  *  M_PI  /  double!  N  ); 
const  double  DEL_Z  =  2.  /  double!  N  ); 

rng:: Random  rng; 

double  x,  y,  phi,  ph2,  z,  zl,  z2,  ph,  th; 
for  (  int  i  =  1;  i  <=  N;  i++  )  { 

z2  =  i  *  DEL_Z; 
zl  =  z2  -  DEL_Z; 

for  (  int  j  =1;  j  <=N;  j++)  { 


ph2  =  j  *  DEL_X; 

phi  =  ph2  -  DEL_X; 

ph  =  rng. uniform!  phi,  ph2  ); 

z  =  rng. uniform!  zl,  z2  )  -  1.; 

//std::cout  «  ph  «  "\t"  «  z  «  std::endl; 
th  =  acos (  z  ) ; 

//std::cout  «  ph  «  "\t"  «  th  «  std::endl; 
x  =  sin(  th  )  *  cos(  ph  ) ; 
y  =  sin(  th  )  *  sin(  ph  ) ; 

std::cout  «  x  «  "\t"  «  y  «  "\t"  «  z  «  std::endl; 

} 

> 

return  EXIT_SUCCESS; 


This  does  a  lot  to  remove  the  clustering  as  shown  in  Fig.  B-3,  but  the  disadvantage  is 
that  we  need  to  know  the  total  number  of  sample  points  beforehand  to  impose  the 
grid.  There  is  another  issue  with  this  stratified  sampling:  to  get  uniform  sampling,  the 
grid  cells  must  be  rectangles  rather  than  squares.  This  results  in  a  different  density 
of  points  along  the  2  dimensions. 

B-3.  Spiral  Distribution 

A  good  discussion  of  the  general  problem  of  distributing  points  uniformly  over 
the  unit  sphere  is  contained  in  the  paper  by  Saff  and  Kuijlaars.3  They  show  that 
a  distribution  of  points  on  the  sphere  spiraling  from  the  north  pole  to  south  pole 
provides  a  good  compromise  that  keeps  the  spacing  between  points  about  the  same. 
Their  formulation  is  implemented  in  Listing  B-3.  This  is  not  randomized,  and  we 
need  to  know  beforehand  the  total  number  of  sample  points.  The  pattern  it  produces 
is  shown  in  Fig.  B-4. 


3Saff  EB,  Kuijlaars  AB.  Distributing  many  points  on  a  sphere.  The  Mathematical  Intelligencer. 
1997;19:5A1. 
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Listing  B-3.  spiral.cpp 


//  spiral.cpp:  uniform  spiral  distribution  over  the  unit  sphere 

//  Implementation  of  the  spiral  distribution  on  the  unit  sphere  as  described  in  the  paper: 

//  Saff  and  Kuijlaars,  "Distributing  Many  Points  on  a  Sphere,"  The  Mathematical  Intelligencer,  Vol.  19  (1967)  pp.  5-11. 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 
using  namespace  std; 

int  main(  int  argc,  char*  argv[]  )  { 

const  double  TW0_PI  =  2.  *  M_PI; 

int  N  =  1000; 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  ); 

double  x,  y,  z,  th,  ph  =  0.; 
for  (  int  i  =  1;  i  <=  N;  i++  )  { 

±f  (  i  —  1  )  { 

z  =  -1.; 
ph  =  0. ; 

} 

else  if  (  i  ==  N  )  { 
z  =  1.; 
ph  =  0. ; 

} 

else  { 

z  =  -1.  +  2.  *  (  i  -  1  )  /  double (  N  -  1  ); 

ph  +=  3.6  /  sqrt(  double (  N  )  *  (  1.  -z)*(l.+z)); 

} 

th  =  acos(  z  ); 

while  (  ph  >  TWO-PI  )  ph  -=  TWO-PI; 

//cout  «  ph  «  "\t"  «  z  «  endl; 

//cout  «  ph  «  "\t"  «  th  «  endl; 

x  =  sin(  th  )  *  cos(  ph  ) ; 

y  =  sin(  th  )  *  sin(  ph  ) ; 

cout  «  x  «  "\t"  «  y  «  "\t"  «  z  «  endl; 

} 

return  EXIT-SUCCESS; 

> 


We  call  this  the  spiral  distribution  since  it  starts  at  the  north  pole  and  spirals  points 
around  the  sphere  until  it  reaches  the  south  pole. 

B-4.  Maximal  Avoidance 


The  last  method  we  consider  is  based  upon  number  theory.  The  code  listing  is  in 
Listing  B-4  and  is  based  upon  the  code  in  Numerical  Recipes .4 


Listing  B-4.  avoidance.cpp 

//  avoidance.cpp:  use  maximal  avoidance  to  generate  a  uniform  distribution  over  the  unit  sphere 

#include  "Random. h" 

#include  "Vector. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

#include  <iomanip> 

va::Vector  spherical!  rng::Random&  rng  )  {  //  returns  a  random  unit  vector  uniformly  distributed  over  the  unit  sphere 

double  x,  y,  z; 

rng . spherical_avoidance(  x,  y,  z  ); 
return  va::Vector(  x,  y,  z  ); 


4Press  WH,  Teukolsky  SA,  Vetterling  WT,  Flannery  BP.  Numerical  recipes  in  C:  the  art  of 
scientific  computing.  New  York  (NY):  Cambridge  University  Press;  1995. 
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int  main(  int  argc,  char*  argv[]  )  { 

unsigned  int  N  =  1024; 

if  (  argc  ==  2  )  N  =  atoi(  argv[l]  ); 

std::cout  «  std : : setprecision (6)  «  std:: fixed; 

rng : : Random  rng; 

for  (  unsigned  int  n  =  0;  n  <  N;  n++  )  std::cout  «  spherical!  rng  )  «  std::endl; 

/* 

double  th,  ph,  xy [2 ] ; 

for  (  unsigned  int  n  =  0;  n  <  N;  n++  )  { 

//rng. avoidance!  xy,  2  ); 

//std::cout  «  xy[0]  *  2.  *  M_PI  «  "\t"  «  xy[l]  *  2.  -  1.  «  std::endl; 
rng. spherical_avoidance(  th,  ph  ); 
std::cout  «  ph  «  "\t"  «  th  «  std::endl; 

} 

*/ 

return  EXIT-SUCCESS; 

} 


This  is  also  deterministic,  not  random,  and  one  requires  to  know  the  number  of 
sample  points  ahead  of  time.  However,  it  does  a  very  nice  job  of  distributing  the 
points,  as  shown  in  Fig.  B-5.  The  points  are  computed  sequentially,  and  number 
theory  is  used  to  avoid  previous  points. 
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Appendix  C.  Some  Properties  of  the  Lognormal  Distribution 
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The  probability  density  function  (PDF)  is 


f(x\fi,a2)  = 


\Z2tt 


exp 


ax 


(In  a;  —  //)" 

2^ 


(C-l) 


and  the  cumulative  distribution  function  (CDF)  is 


F(x  \/x,a2)  =  -[1  +  erf 


in  x  —  /i 

V2a 


(C-2) 


where  /i  is  the  location  parameter  and  a  is  the  scale  parameter.  Expressions  for  the 
usual  metrics  are  listed  in  Table  C-l. 


Table  C-l.  Properties  of  the  lognormal  distribution 


Statistic 

Expression 

Geometric  Mean 

xg  = 

Geometric  Standard  Deviation 

<7g  =  ea 

Median 

x50  = 

Mean 

x  = 

Mode 

x  =  e^-*2 

•  The  values  ax  and  x / a  are  equally  likely  for  any  value  a  ^  0.  That  is, 

f(ax)  =  f(x/a). 

•  68%  of  the  distribution  is  contained  in  the  interval  [xgag  xgag\. 

•  95%  of  the  distribution  is  contained  in  the  interval  [xgag  2,  xga2]. 


The  n-th  moment  about  the  origin  is 

/»oo 

An  =  /  xnf(x)dx 


\[2i ra  Jo 

i  r 


enlnxf(x)dx 

I'OO 

exp 

exp 

=  exp  (  n/j  +  -nV2 


=  x'a  exp  I  ~n2 a2 


n a  Jo 


( 1  TO  T*  -  /  /  ^  2 

1  -  ,/j  +  n  In  x 
2a- 

(In  x  —  n  —  na2)2 


d:  In  x 


2a 2 


+  nix  +  -n2<r2 


d  In  x 


(C-3) 
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The  n-th  moment  distribution  function  of  f(x  |  /i,  o'2)  is  defined  by 


fn(x\n,a2) 


—xnf(x\n,a2). 

An 


(C-4) 


Using  Eqs.  C-3  and  C-l, 

t  I  ..  _2\ _ _ 


1  2 

I  1 

n(i  — 

-n  (J 

2  2 

—  exp 
\2ttctx 

1  2  2) 

1  1 

nil  — 

-n  a 

2  2 

/ST^  eXP 
\J2/nax 

-  exp 

(In  2; 

—  [1  —  TUT2)2 

2  a2 

(In  a;  —  nY 


+  n  In  x 


2a2 

(In  a;  —  fi  —  na2)2 
2^2 


+  n/a  +  -n2a2 


(C-5) 


And  thus  we  have  derived  the  Fundamental  Theorem  of  the  Moment  Distribution'. 


The  n-th  moment  distribution  of  a  lognormal  distribution  with  parameters  //  and  a2 
is  also  a  lognormal  distribution  with  parameters  /i  +  no2  and  a2,  respectively, 

fn(x  |  /x,  a2)  =  f(x  |  n  +  no-2,  o'2).  (C-6) 


We  can  also  show  that  the  product  and  quotient  of  lognormal  distributions  are  also 
lognormal.  Using  the  notation  of  Aitchison  and  Brown1,  if  X1  ~  A(//| .  of)  and 
X2  ~  A (/i2,  erf),  then  the  product  X\X2  is  also  lognormal  with 


XiX2  ~  A(/i!  +  fi2 ,  o’2  +  o'2), 


(C-7) 


and  the  quotient  X\/ X2  is  also  lognormal  with 

X\jX2  ~  A(/i!  —  /i2, cr2  +  0-2). 


(C-8) 


We  can  derive  these  results  as  follows. 


'Aitchison  J,  Brown  JAC.  The  lognormal  distribution.  New  York  (NY):  Cambridge  University 
Press;  1963. 


Approved  for  public  release;  distribution  is  unlimited. 


79 


C-1.  Product  of  Lognormals 


Let  X\  ~  A(//| .  erf)  and  X2  ~  A(/i2,<t|)  be  lognormal  distributions.  Then  the 
cumulative  distribution  of  their  product  is 


F 


Xi_X2 


u  = 


fix  i,  x2)  dx  i  dx  2 


Xi.Y2<« 


f(xi,x2)  dx 


dx  i- 


(C-9) 


The  density  is  obtained  by  differentiating  with  respect  to  u,  so  that 


fx1x2(u)  =  [  f  (xi,—')  —dxi 
Jo  V  X\ )  X\ 

f°°  Ju\ 

=  /  f(x)f  (  —  )  din  a; 

Jo  ' x ' 

(lna;-^i)2l  1 


r°°  i 

/ o  s/Fk  o\X 
1 


exp 


J  ^/2tt a2(u/x) 


exp 


(ln(w/a;)  -  fi2) 
2  a\ 


21 


din  a; 


2nai<j2u  Jo 

-J—f 

2TTa1a2u  Jo 

1  r 


2ira1a2u  J_c 

i  r 


exp 


exp 


exp 


(In  a;  —  /zi)2  (In  u  —  lna;  —  n2f 


2<rf 


2  cr| 


din  a: 


(lna:  — /.ii)2  (lna:  +  n2  —  In  it)2 
2erf  2ct| 


din; 


(x  —  /ii)2  (x  +  /j.2  —  lnw)2 


27r<Ticr2'u  J_c 


2crf 

exp  [-ff(x)]  dx, 


2ct| 


da; 


(C-10) 


where 

anti 


-  {x  +  XlXixFX  ^  a(x  _  ^  _  c)=,  (c-ii) 


2<rf 


2a| 


a=^,  b  =  c=lnu. 


2  erf 7  2  erf 


(C-1 2) 


Completing  the  square  in  x  gives,  after  some  messy  algebra, 


g(x)  =  (a  +  b)  (  x  -  — — b^2  +  bC^j  + -^y  [c  -  (;Ui  + /i2)]2  (C-1 3) 

a  +  b  J  a  +  b 


Approved  for  public  release;  distribution  is  unlimited. 


80 


or,  in  terms  of  0\ ,  <r2,  and  In  u, 


„(„\  _  ai  +  °2  („  a2Hi  ~  aid2  +  of  lnufz  |  (lnw  -  Hi  —  /a2)s 
fix3')  —  0  2  2  I  x  o-  o  I  + 


2(Tjcr2 

Therefore,  returning  to  Eq.  C-10, 

(Inw  -  AO  -  M2)2 


af  +  al 


2K  +  af) 


1x^2  ( U )  =  exp 


1 


27T(Tl(J2M  7- 


2(cri  +  o'2) 
1*00 

exp 


erf  +  CTo  /  CT9M1  —  ctiM2  +  erf  lnu 

X - - - o - o - 


2crjcr. 


2^2 
2 


err  +  crS 


The  integral  is  now  easily  evaluated  with  the  substitution 


£  = 


Vai  + 1 
72 


erier2 


&2H1  —  +  cr2  lnw 


a'l  +cr| 


and  d£  = 


7°~2  + 1 

72 


da; 


ericr2 


and  we  get 

fx1x2(u)  =  exp 


(In  u  —  Hi  -  /i2)s 

2K+^f) 


,-7 


(C-14) 


dx  (C- 15) 


(C- 16) 


a/27t(o-2  +  of)  u 


exp 


72  'K\Jc>\  +  of  U  J-C 
(Inw  -  Hi  —  V2)2 

2KTo|)  _ 


— A(/o  +  H2,  of  +  of). 


(C-17) 


Thus, 

X1X2  ~  A(/ri  +  H2-,  of  +  of) 


as  was  to  be  shown. 


(C-18) 


C-2.  Quotient  of  Lognormals 

Let  X\  ~  A(//| .  er2)  and  AT  ~  A(/e2,  of)  be  lognormal  distributions.  Then  the 
cumulative  distribution  of  their  quotient  is 


FXl/x2(u)  = 


f(x  1,  x2)  dx\  dx2 


X1/X2<u 


f(x i,x2)  dx 


dx2- 


(C-19) 
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The  density  is  obtained  by  differentiating  with  respect  to  u,  so  that 


nOO 

fx  1  /x2  (u)=  f  ( ux2 ,x2)x2  dx2 

Jo 

/oo 

f(ux)f(x)  x  dx 


V2n 


G\X 


exp 


(ln(ux)  -  Hi)2 


■  exp 


(In a:  -  H2)2 


1 


2ir<Jia2u  J0 

1  r 


2ir<Jia2u  J0 

i  r 


2nai<j2u  J_ c 

i  r 


exp 


exp 


exp 


2c2  \  y/2TT  cr2x 

(In  it  +  lnx  —  Hi)2  (lnx  — /X2)2"1 


2  a\ 


2cr| 


2a\ 


d  In: 


(lnx  —  Hi  +  In  it)2  (lnx  —  ^t2)2x  1 


2  cr2 


2a2 


d  In  x 


(x  —  hi  +  lntt)2  +  (x  —  H2)2 


2a'l 


2<t\ 


dx 


2nai<j2u  J_c 


exp  [— h(x)\  dx, 


where 
h(x) 
and 


_(x-fl1  +  \n  u)2  (x-/i2)2_  V2  ,  ,/  n2 

- — ^ - t - — 2 -  =  0(3?  —  /Xi  +  c)  +  6(x  —  H2)  1 

Act  -y  2j(J  <2 


a=2^’  C"ln“- 


Completing  the  square  in  x  gives 


h(x)  =  (a  +  b)  (  x 


d./H  T  b/ji 2  —  oc 


a  +  b 


+  _^|c-(,1-,2)|2 


or,  in  terms  of  <Ji,  cr2,  and  In  u, 


,  ,  ,  cr2  +  cr| 

^'(•r)  —  o_2^2 


2cr{cr. 


cr|/ri  +  a2;U2  -  cr|  In  u  \  [In  u  -  (/ii  -  A^)]5 
X  o  ;  q  I  T 


<^1+^2 


2(cri  +  ^22) 
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x  dx 

(C-20) 

(C-21) 

(C-22) 

(C-23) 

(C-24) 


Therefore,  returning  to  Eq.  C-20, 


fx2/x1(u)  =exp 


[In  it  -  {hi  - 


1 


27T£ricr2tt  J_ 


2(al+al) 

»oo 

exp 


a i  +  &2  f  +  tTi/i2  —  cr2  In  w 


2^.2 
2 


2(7^  (7 . 


(7i  +cr| 


dx  (C-25) 


The  integral  is  now  easily  evaluated  with  the  substitution 


£  = 


\/crj  +  ^  cr|/ii  +  (7i/x2  —  tr|  lnu 

\/2crit72 


a'l  +a| 


and  d£  = 


Vai  + 1 
72 


CTit72 


dx  (C-26) 


and  we  get 

/x2/x,  (w)  =  exp 


[In  w  -  (n i  -  /x2)]: 


e 


2(^1  +  cr|)  J  v^vta/o-i  +  of 


u  J- c 


a/27t((T^  +  <r| 


u 


exp 


[In  w  -  (/ii  -  /x2)p 

2(^  +  *22) 


=A(/il  -fJL2,Oi  +  o|). 


(C-27) 


Thus, 

X1/X2  ~  A(/xi  —  /r2,  of  +  o^) 


as  was  to  be  shown. 


(C-28) 


In  general,  if  Xt,  i  =  1,2,  ...,n  are  independent  random  variables  with  Xt  ~ 
A(/Xj,  a,?)  and  a  and  p,  are  constants,  then2 


n  n  n 

a  n  Xi  i  ~  A(ln  a  +  55  51  pia^ 

i=l  i= 1  2—1 


(C-29) 


Notice  that  the  vai'iance  can  only  increase,  never  decrease.  The  next  section  gives  an 
application  of  this  result. 


2Crow  EL,  Shimizu  K,  Editors.  Lognormal  distributions:  theory  and  applications.  New  York 
(NY):  Marcel  Dekker,  Inc.;  1988. 
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C-3.  Mass  per  Unit  Area  Distribution 


Suppose  that  the  dimensionless  shape  factor  7  is  lognormally  distributed  with  7  ~ 
A  (hi,  a2)  and  that  fragment  mass  is  also  lognormally  distributed  with  m  ~  A  (/r2,  07). 
Then,  since  the  fragment  presented  area,  Ap,  equals  7 (m/p)2/3,  where  p  is  the 
material  density,  it  follows  from  Eq.  C-29  that  mass  per  unit  area  is  also  lognormally 
distributed: 

^-  =  V/3m1/3~  Afea2),  (C-30) 

Ap  7 

where 

P  =  hi  p2/3  -  pi  +  y  and  a2  =  a2  +  (y)  .  (C-31) 

To  illustrate  this,  let  us  return  to  the  artillery  fragments  described  in  Section  3.2. 
A  lognormal  distribution  fits  the  fragment  masses  (in  grams)  with  //2  =  1.690  and 
<7 2  =  1.323.  We  also  found  that  the  shape  factor  distribution  was  lognormal  with 
Pi  =  0.597  and  a  1  =  0.341.  From  Eqs.  C-29,  C-30,  and  C-31,  it  then  follows  that 
the  mass  per  unit  area  distribution  (in  g/cm2)  should  also  be  lognormal  with 

p  =  In  7.8 32/3-p1+p2/3  =  1.338  and  a  =  ^  a2  +  (a2/3)2  =  0.557.  (C-32) 

The  program  in  Listing  C-l  generates  the  shape  factor  and  mass  independently  and 
outputs  the  mass  per  unit  area. 


Listing  C-l.  rnu.cpp 

//  mu.cpp:  generate  shape  factor  and  mass  independently 

//  and  output  mass  per  unit  area  to  compare  to  the  theoretical  distribution 

#include  "Random. h" 

#include  <iostream> 

#include  <cstdlib> 

#include  <cmath> 

int  main(  void  )  { 

const  int  N  =  10000; 

const  double  MU1  =  0.597,  SIGMAl  =  0.341,  RH0  =  7.83,  C  =  pow(  RH0,  2./3.  ); 
const  double  MU2  =  1.690,  SIGMA2  =  1.323; 
rng : : Random  rng; 
double  sf,  m,  mu; 

for  (  int  i  =  0;  i  <  N;  i++  )  { 

sf  =  rng . lognormal (  0.,  MU1,  SIGMA1  ); 
m  =  rng.lognormal(  0.,  MU2,  SIGMA2  ); 
mu  =  C  *  pow(  m,  1./3.  )  /  sf; 
std::cout  «  mu  «  std::endl; 

> 

return  EXIT_SUCCESS; 


//  shape  factor  (dimensionless) 
//  mass  (grams) 

//  mass  per  unit  area  (g/cnT2) 


These  values  are  compared  to  the  theoretical  distribution  in  Fig.  C-l. 
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Fig.  C-l.  The  program  in  Listing  C-l  was  used  to  generate  independent  samples  of  shape  factor 
and  mass  and  output  mass  per  unit  area  in  order  to  compare  to  the  theoretical  distribution 
(black  curve) 


This  shows  that  the  resulting  distribution  is  indeed  lognormal  with  fi  and  a  as  given 
by  Eq.  C-3 1 .  Thus,  the  mass  per  unit  area  samples  could  be  generated  from  a  single 
lognormal  distribution. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 

TERMS: 

3D:  3  dimensional 

FATEPEN:  Fast  Air  Target  Encounter  Penetration 
RCC:  right-circular  cylinder 

RPP:  rectangular  parallelepiped  (also  known  as  cuboid) 
STL:  stereolithography 
MATHEMATICAL  SYMBOLS: 

/,  PDF:  probability  density  function 
F,  CDF:  cumulative  distribution  function 
7:  shape  factor 
p:  material  density 

A(/z,  a2):  lognormal  distribution  with  mean  //  and  variance  a2 
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